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1.  Introduction 

These  notes  are  arranged  in  the  following  manner.  In  the  introduction  section,  we 
use  several  examples  to  illustrate  some  of  the  issues  that  must  be  addressed  when  we 
model  hypersonic  flows.  This  leads  to  a  discussion  of  what  types  of  computational  fluid 
dynamics  methods  are  suitable  for  these  flows.  Then  the  conservation  equations  for  a 
mixture  of  chemically  reacting  and  weakly  ionized  gases  is  developed.  We  discuss  the 
thermochemistry  models  and  the  relevant  boundary  conditions  for  these  flows.  Then  in 
the  third  section,  computational  fluid  dynamics  methods  for  these  flows  are  discussed.  We 
analyze  the  conservation  equations,  and  discuss  an  upwind  method.  Then,  the  integration 
of  the  source  terms  is  discussed.  In  the  fourth  section  we  discuss  several  advanced  topics 
in  the  modeling  of  hypersonic  flows. 

1.1  Examples  of  Hypersonic  Flows 

In  this  subsection,  we  discuss  a  few  external  hypersonic  flows  to  illustrate  some  of  the 
issues  that  we  must  address  when  the  flow  field  is  chemically  reacting.  These  examples 
show  that  aerothermochemistry  can  have  a  major  impact  on  aerodynamic  coefficients,  heat 
transfer  rates,  and  radiative  emission  from  hypersonic  flows. 

1.1.1  Hypersonic  Double-Cone  Flow 

An  interesting  flow  field  is  created  by  a  double-cone  geometry.  Figure  1.1  plots  a 
schematic  of  this  flow  field  for  an  approximately  Mach  12  free-stream  condition.  Note  the 
attached  shock  wave  that  originates  at  the  first  cone  tip,  the  detached  shock  wave  formed 
by  the  second  cone,  and  the  resulting  shock  triple  point.  The  transmitted  shock  impinges 
on  the  second  cone  surface,  which  separates  the  flow  and  produces  a  large  localized  increase 
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in  the  pressure  and  heat  transfer  rate.  This  pressure  rise  causes  the  flow  to  separate,  and 
also  produces  a  supersonic  under-expanded  jet  that  flows  downstream  near  the  second  cone 
surface.  The  size  of  the  separation  zone  depends  strongly  on  the  location  and  strength  of 
the  shock  impingement.  This  flow  field  is  very  sensitive  to  the  wind  tunnel  conditions,  the 
physical  models  used  in  the  CFD  code,  and  the  quality  of  the  numerical  methods  used  to 
predict  the  flow  (Ref.  40). 


Fig.  1.1  -  Schematic  of  the  hypersonic  double-cone  flow  field. 

Experiments  on  the  double-cone  have  been  performed  at  hypersonic  conditions  in 
the  CUBRC  Large  Energy  National  Shock  Tunnel  (LENS).  These  experiments  used  a 
large  model  with  many  surface-mounted  heat  flux  and  pressure  transducers  (Ref.  22). 
Nitrogen  was  used  as  the  test  gas  to  minimize  the  effects  of  chemical  reactions,  and  the 
experiments  were  done  at  low  density  to  ensure  laminar  boundary  layers  and  shear  layers. 
This  dataset  was  the  subject  of  a  blind  code  validation  study.  In  general,  the  comparison 
between  simulation  and  experiment  was  shown  to  be  quite  good,  however  there  were  several 
important  differences.  Interestingly,  the  simulations  performed  with  high-quality  numerical 
methods  on  the  finest  grids  slightly  over-predicted  the  size  of  the  separation  zone,  and  all 
simulations  predicted  excessive  heating  in  the  attached  region  prior  to  separation  (Ref.  22). 
This  is  shown  in  Fig.  1.2,  which  presents  typical  results  for  two  double-cone  cases.  The 
error  on  the  first  cone  is  as  much  as  20%,  which  is  particularly  puzzling  because  the 
pressure  is  accurately  predicted  in  this  region.  Many  attempts  were  made  to  explain  this 
difference  by  running  CFD  cases  with  extreme  grid  resolution,  finite  nose-tip  bluntness, 
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model  misalignment,  and  uncertainties  in  reaction  rates.  None  of  these  effects  was  found 
to  explain  the  differences  shown  in  Fig.  1.2. 


Fig.  1.2  -  Comparison  of  predicted  and  measured  surface  pressure  and  heat  flux  on  the 
double  cone  model.  Ref.  22. 


Fig.  1.3  -  Comparison  of  measured  double  cone  heat  flux  and  different  models  for  the 
effects  of  nonequilibrium  vibrational  energy.  ’Nominal’  is  the  baseline  model  (Fig.  1.2); 
’Noneq.  ’  is  accounting  for  vibrational  nonequilibrium  in  the  free-stream;  and  ’Noneq.  Slip  ’ 
also  includes  the  effect  of  imperfect  accommodation  of  vibrational  energy  to  the  surface. 
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The  specification  of  the  free-stream  conditions  in  a  hypersonic  shock  tunnel  can  be 
difficult  because  these  facilities  may  be  subject  to  non-ideal  effects  in  the  nozzle.  Namely, 
a  reflected  shock  wave  is  used  to  heat  and  compress  the  test  gas  to  extremely  high  pressure 
and  temperature,  which  results  in  vibrational  excitation  and  chemical  reaction.  Then  the 
test  gas  rapidly  expands  through  the  nozzle,  and  its  internal  energy  state  may  not  fully 
de-excite  during  the  expansion.  As  a  result,  the  gas  flowing  over  the  model  may  be  in 
a  non-ideal  thermo-chemical  state.  The  double-cone  experiments  were  run  with  nitrogen 
at  moderate  enthalpy  (Jiq  <  4MJ/kg);  this  results  in  virtually  no  chemical  reaction,  but 
vibrational  excitation  of  the  gas  in  the  reflected  shock  region.  Nitrogen  vibrational  modes 
relax  very  slowly,  and  for  these  test  conditions  this  results  in  elevated  vibrational  energy 
in  the  wind-tunnel  test  section.  A  vibrational  finite-rate  simulation  of  the  nozzle  flow 
shows  that  the  vibrational  energy  modes  are  frozen  near  the  throat  temperature  (Tv  — 
2560  K).  This  has  two  major  effects:  the  kinetic  energy  flux  is  reduced  by  about  10%, 
and  because  nitrogen  vibrational  energy  modes  are  inefficient  at  accommodating  to  most 
metallic  surfaces,  they  do  not  transfer  their  energy  to  the  model.  These  two  effects  reduce 
the  heat  flux  by  about  20%,  and  significantly  improve  the  comparison  between  CFD  and 
experiment;  this  is  shown  in  Fig.  1.3  (Ref.  40). 

1.1.2  Mars  Science  Laboratory  Entry 

During  the  past  20  years  or  so,  great  strides  have  been  made  in  the  simulation  of 
atmospheric  entry  flows.  In  part,  this  is  a  result  of  increases  in  computer  power  and  the 
development  of  efficient  parallel  codes  for  solving  the  governing  equations.  Recently,  un¬ 
structured  grids  have  begun  to  be  used  for  difficult  aerothermodynamics  problems.  Inter¬ 
estingly,  however  almost  no  research  has  been  done  on  improving  the  governing  equations 
used  in  these  simulations.  Thus,  we  are  solving  essentially  the  same  equations  now  as  when 
the  first  CFD  simulations  of  nonequilibrium  atmospheric  entry  flows  were  performed.  Fur¬ 
thermore,  current  codes  are  often  used  far  beyond  the  range  of  conditions  for  which  they 
have  been  validated.1 

Figure  1.4  illustrates  a  recent  simulation  of  the  flow  over  the  Mars  Science  Laboratory 
(MSL)  capsule  at  Mach  18.1  (Ref.  53),  with  the  flow  in  the  wake  region  represented  by 
the  detached  eddy  simulation  (DES)  approach  (Ref.  59)  with  the  Spalart-Allmaras  one- 
equation  turbulence  model  (Ref.  58)  and  the  Catris  and  Aupoix  (Ref.  7)  compressibility 
correction.  This  simulation  used  finite-rate  kinetics  to  represent  the  Mars  atmosphere 

1  The  issue  of  code  validation  is  beyond  the  scope  of  these  notes,  but  it  is  particularly  difficult  for 
hypersonic  flows  because  it  is  very  difficult  to  isolate  specific  elements  of  the  physical  models  in 
well-controlled  experiments. 
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and  was  performed  on  16  million  element  hybrid  unstructured  grid.  Not  shown  in  this 
image  is  that  the  operation  of  the  reaction  control  system  thrusters  is  also  modeled  in  this 
simulation.  With  current  parallel  clusters,  such  a  calculation  can  be  performed  at  useful 
turn-around  times.  (The  base  flow  is  established  in  about  half  a  day,  and  the  operation  of 
the  RCS  thrusters  is  computed  and  averaged  over  an  additional  two  days  on  240  cores  of 
an  SGI  Altix  SE  1300  cluster.) 


Fig.  1.4  -  Temperature  contours  in  the  flow  field  of  the  MSL  capsule  at  Mach  18.1  condi¬ 
tions. 


1.2  Important  Effects 

1.2.1  Thermochemical  Nonequilibrium 

A  gas  is  in  thermal  nonequilibrium  if  its  internal  energy  cannot  be  characterized  by 
a  single  temperature,  and  it  is  in  chemical  nonequilibrium  if  its  chemical  state  does  not 
satisfy  chemical  equilibrium  conditions.  Portions  of  many  external  hypersonic  flows  are  in 
thermal  and  chemical  nonequilibrium.  This  occurs  because  as  the  gas  passes  through  the 
bow  shock  wave,  much  of  its  kinetic  energy  is  converted  to  random  translational  motion. 
Then,  collisions  transfer  translational  energy  to  rotational,  vibrational,  electronic,  and 


RTO-EN-AVT-162 


15-5 


Computational  Fluid  Dynamics  for  Atmospheric  Entry 


ORGANIZATION 


chemical  energy.  This  energy  transfer  takes  a  certain  number  of  collisions,  during  which 
time  the  gas  moves  to  a  new  location  where  the  temperature  and  density  may  be  different. 
Thus,  the  internal  energy  modes  and  chemical  composition  of  the  gas  lag  the  changes  in 
the  translational  temperature.  We  can  determine  if  a  flow  will  be  in  thermal  or  chemical 
nonequilibrium  by  constructing  the  Dahmkohler  number,  Da,  which  is  the  ratio  of  the 
fluid  motion  time  scale  to  the  internal  energy  relaxation  or  chemical  reaction  time  scale. 

Consider  the  steady-state  mass  conservation  equation  for  species  s 


d 

dxj 


(psUj) 


=  Wa, 


(1.1) 


where  ps  is  the  mass  density  of  species  s,  Uj  is  the  gas  velocity  in  the  Xj  direction,  and  ws 
is  the  rate  of  production  of  species  s  per  unit  volume  due  to  chemical  reactions.  We  can 
non-dimensionalize  this  equation  using  the  total  density,  p,  the  speed,  V,  and  a  relevant 
length  scale  such  as  the  nose  radius,  rn.  Then 


d  ,  _  x  rnws 

DD  -  pV 


(1.2) 


Thus,  Da  represents  the  ratio  of  the  fluid  motion  time  scale,  Tf,  to  the  chemical  reaction 
time  scale,  rc;  or  it  is  the  ratio  of  the  chemical  reaction  rate  to  the  fluid  motion  rate.  A 
similar  expression  can  be  derived  to  describe  the  relative  rate  of  internal  energy  relaxation. 

When  Da  — >  oo,  the  internal  energy  relaxation  or  chemical  reaction  time  scale  ap¬ 
proaches  zero  (becomes  infinitely  fast),  and  the  gas  is  in  equilibrium.  That  is,  its  thermal 
or  chemical  state  adjusts  instantaneously  to  changes  in  the  flow.  When  Da  — >  0,  the  reac¬ 
tion  time  scale  approaches  infinity,  the  gas  is  frozen  and  does  not  adjust  to  changes  in  the 
flow.  The  Dahmkohler  number  is  useful  for  determining  how  reactive  the  gas  is,  and  what 
type  of  analysis  is  appropriate  for  given  flow  conditions. 

When  the  chemical  source  term,  vus,  is  proportional  to  the  density  squared  (as  it  is  for 
dissociation  reactions),  the  binary  scaling  law  can  be  derived  from  the  above  expression. 
Let  us  write  ws  —  Cp2kf,  where  kf  is  a  temperature-dependent  reaction  rate,  and  C  is  a 
constant.  Then  we  have 

Da  —  prn^^~.  (1.3) 

Thus,  the  reaction  rate  is  proportional  to  the  density- length  scale  product,  kf  depends 
exponentially  on  temperature,  which  for  hypersonic  flows  depends  on  the  free-stream  ki¬ 
netic  energy  Therefore,  for  a  dissociation-dominated  flow,  the  Dahmkohler  number 

depends  on  the  binary  scaling  parameter,  prn,  and  the  free-stream  kinetic  energy. 
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1.2.2  Vibration-Dissociation  Coupling 

When  a  gas  becomes  vibrationally  excited,  the  population  of  the  excited  vibrational 
states  increases.  As  shown  in  Fig.  1.5,  this  decreases  the  energy  required  to  dissociate  the 
molecule.  Therefore,  the  vibrational  state  of  a  molecule  affects  its  dissociation  rate.  This 
process  is  not  fully  understood,  and  simple  models  that  can  be  implemented  in  computa¬ 
tional  methods  are  largely  unvalidated.  Even  small  changes  in  the  dissociation  rate  can 
change  the  flow  field  considerably  and  can  lead  to  uncertainties  in  the  trim  angle  of  attack 
of  hypersonic  vehicles. 

There  are  many  models  available  for  the  vibration-dissociation  process,  and  we  will 
discuss  the  most  popular  models  in  these  notes. 


Fig.  1.5  -  Schematic  of  the  vibrational  state  of  a  molecule  and  the  energy  required  for 
dissociation. 


1.2.3  Finite-Rate  Wall  Catalysis 

One  of  the  most  important  parameters  that  determines  the  convective  heat  transfer 
rate  for  hypersonic  vehicles  is  the  surface  catalytic  efficiency.  Fay  and  Riddell  (Ref.  11) 
used  a  self-similar  stagnation  point  boundary  layer  analysis  to  show  that  depending  on  the 
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reactivity  (Dahmkohler  number)  of  the  boundary  layer  and  the  catalytic  efficiency  of  the 
surface,  the  heat  transfer  rate  is  dramatically  changed.  As  seen  in  Fig.  1.6,  if  the  boundary 
layer  is  frozen  (Da  — >  0)  and  the  body  surface  is  noncatalytic,  the  heat  transfer  rate  may 
be  reduced  by  50%  or  more,  depending  on  the  fraction  of  energy  tied  up  in  the  chemical 
energy  modes. 


Fig.  1.6  -  Effect  of  wall  catalysis  on  convective  heat  transfer  rate  at  a  stagnation  point. 
( After  Fay  and  Riddell  (Ref.  11).) 

When  the  wall  is  noncatalytic,  it  does  not  promote  recombination  at  the  surface.  Thus, 
if  the  reaction  rate  is  slow  near  the  surface,  a  fraction  of  the  total  energy  of  the  gas  remains 
in  the  form  of  chemical  energy.  Thus,  it  does  not  contribute  the  convective  heating.  If, 
on  the  other  hand,  the  surface  is  catalytic  or  the  boundary  layer  is  near  equilibrium,  the 
chemical  energy  is  released  at  the  vehicle  surface,  and  the  heat  transfer  increases. 

It  is  important  to  note  how  the  chemical  reactions  scale  with  density.  As  we  saw  in 
the  previous  section,  dissociation  scales  with  the  binary  scaling  parameter,  prn.  However, 
recombination  is  a  three-body  process,  and  as  a  result  scales  with  p2 .  Thus,  in  typical 
low-density  atmospheric  entry  conditions,  recombination  may  be  very  slow  relative  to 
dissociation.  The  surface  reactions  typically  scale  differently  because  they  are  usually 
diffusion  limited  and  depend  strongly  on  the  surface  material  properties. 

1.2.4  Nonequilibrium  Thermal  Radiation 

The  modeling  of  thermal  radiation  from  the  flow  field  remains  a  major  challenge.  The 
difficulty  arises  because  the  radiative  emission  from  the  gas  depends  very  strongly  on  its 
internal  energy.  For  example,  immediately  behind  the  bow  shock  wave  in  the  stagnation 
region  of  a  vehicle,  the  vibrational  temperature  may  overshoot  the  equilibrium  post-shock 
temperature.  Then,  if  the  population  of  the  excited  electronic  states  of  the  gas  is  governed 
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by  the  vibrational  temperature,  there  is  a  dramatic  super-equilibrium  population  of  the 
excited  electronic  states.  The  excited  electronic  states  decay  to  the  ground  state,  and  the 
electronic  energy  is  emitted  as  photons  (thermal  radiation),  many  of  which  are  absorbed  by 
the  body  surface.  Thus,  for  high-energy  flows  where  radiative  heating  is  important,  there 
may  be  a  significant  increase  in  the  heat  transfer  rate  due  to  nonequilibrium.  This  process 
is  very  difficult  to  model  because  there  are  many  complicated  rate-dependent  processes 
competing  for  the  thermal  energy  produced  by  the  shock  wave. 

1.2.5  Low  Density  Effects 

During  re-entry  many  hypersonic  vehicles  spend  a  significant  time  at  high  altitude 
where  the  gas  dynamics  are  poorly  understood.  In  this  high  Knudsen  number  regime,  the 
bow  shock  wave  thickness  becomes  an  appreciable  fraction  of  the  shock  standoff  distance, 
none  of  the  internal  energy  modes  are  in  equilibrium  with  the  translational  modes,  and 
there  is  velocity  and  temperature  slip  at  the  body  surface.  Also,  the  chemical  reactions 
may  not  conform  to  standard  reaction  rate  models  and  the  vibration-dissociation  coupling 
effects  are  very  important.  Reasonable  progress  has  been  made  in  the  modeling  of  each 
these  effects  (Refs.  5,  14,  32).  However,  in  extremely  low  density  flows,  the  continuum 
formulation  completely  breaks  down  and  a  particle  based  simulation  method  must  be  used 
(such  as  the  direct  simulation  Monte  Carlo  method  (Ref.  2)). 

1.2.6  Other  Effects 

Many  other  complicated  phenomena  may  occur  in  hypersonic  external  flows.  For  ex¬ 
ample,  at  low-earth-orbit  re-entry  speeds  and  above,  the  flow  field  becomes  ionized  and 
the  radio  frequency  transmissions  may  be  blacked  out.  At  higher  speeds,  the  vehicle  sur¬ 
face  must  be  ablative  to  protect  the  vehicle  from  heating.  In  this  case,  foreign  species 
are  injected  into  the  flow  field,  where  they  react  with  the  air  species  to  form  other  chem¬ 
ical  species.  Thirdly,  shock-wave  boundary  layer  interactions  are  very  complicated  and 
intense  at  hypersonic  conditions.  Finally,  of  course,  there  is  the  question  of  transition  to 
turbulence.  If  the  boundary  layer  is  turbulent,  the  convective  heat  transfer  rate  to  small- 
angle  bodies  (blunted  cones,  for  example)  increases  by  a  factor  of  between  three  and  eight. 
There  are  no  reliable  models  for  transition  at  hypersonic  Mach  numbers,  and  the  effect  of 
boundary  layer  chemical  reactions  on  transition  is  poorly  understood.  However,  it  has  been 
shown  experimentally  and  with  stability  theory  that  endothermic  (energy  absorbing)  reac¬ 
tions  tend  to  delay  transition  in  hypersonic  boundary  layers  (Ref.  23).  Reynolds-averaged 
Navier-Stokes  (RANS)  models  are  largely  unproven  at  hypersonic  conditions. 
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1.3  Computational  Methods 

From  the  previous  discussion,  it  is  clear  that  there  are  many  difficult  issues  that  must 
be  addressed.  Hypersonic  flows  may  have  one  or  more  of  the  following  features: 

1.  A  large  number  of  chemical  species  reacting  with  one  another  at  a  wide  range  of  relative 
time  scales. 

2.  Complicated  interactions  between  the  fluid  motion,  internal  energy  state,  and  the 
chemical  composition  of  the  gas. 

3.  Length  scales  ranging  from  the  characteristic  length  of  the  body  down  to  the  shock 
wave  thickness. 

4.  Complex  gas-surface  interactions,  including  slip  at  the  wall,  foreign  species  injection, 
and  finite-rate  catalysis  of  reactions  at  the  surface. 

5.  Transitional  and  turbulent  boundary  layers. 

A  numerical  method  that  can  solve  the  equations  that  describe  these  flows  must  have 
certain  qualities.  Generally,  we  are  interested  in  steady-state  flows.  Also,  many  flows  have 
regions  of  very  high  reaction  rate,  making  the  range  of  time  and  length  scales  very  large. 
Thus,  it  is  essential  to  use  an  implicit  method  for  at  least  the  chemical  reaction  terms. 
Also,  the  computational  and  memory  costs  of  the  method  should  not  increase  too  quickly 
with  the  number  of  chemical  species  being  considered.  Finally,  to  solve  very  large  and 
difficult  problems,  it  is  mandatory  that  the  method  run  efficiently  on  parallel  computers. 
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2.  Conservation  Equations 


2.1  Assumptions 

We  assume  that  the  gas  is  described  by  the  Navier-Stokes  equations  extended  to 
account  for  the  presence  of  chemical  reactions  and  internal  energy  relaxation.  For  these 
equations  to  be  valid,  the  flow  must  satisfy  the  following  criteria: 

1.  The  gas  must  be  a  continuum.  If  we  relate  the  mean-free-path,  A,  to  a  local  length  scale 
that  is  determined  by  the  normalized  density  gradient,  we  can  form  the  gradient-length- 
local  Knudsen  number  as  (Ref.  4): 


(Kn)GLL 


A  dp 
p  di 


(2.1) 


where  t  is  in  the  direction  of  the  steepest  density  gradient.  Boyd,  et.  al  (Ref.  4)  showed 
that  when  (Kii)Gll  >  0.05  the  Navier-Stokes  formulation  fails.2  This  typically  occurs  in 
shock  waves  and  near  the  body  surface  in  low  density  flows.  Also,  the  wake  region  of 
blunt  bodies  may  have  regions  of  continuum  formulation  failure.  In  these  regions,  either  a 
higher-order  continuum  formulation  must  be  used  (Refs.  32,  71)  or  a  particle-based  method 
such  as  the  direct  simulation  Monte  Carlo  (DSMC)  method  is  required. 

2.  The  mass  diffusion  fluxes,  shear  stresses,  and  heat  fluxes  must  be  proportional  to  the 
first  derivatives  of  the  flow  quantities.  If  not,  a  non-continuum  approach,  such  as  the  direct 
simulation  Monte  Carlo  (DSMC)  method  must  be  used.3 

3.  The  internal  energy  modes  must  be  separable.  That  is,  we  can  describe  each  energy 
mode  by  its  own  temperature.  For  example,  the  energy  contained  within  the  vibrational 
modes  cannot  be  a  function  of  the  rotational  energy  state. 

4.  Finally,  the  flow  is  only  weakly  ionized.  In  this  case,  the  Coulomb  cross-section  is  small 
relative  to  the  electron-neutral  collision  cross-section. 


2  In  this  work,  failure  was  defined  to  occur  when  the  Navier-Stokes  solution  differs  by  5%  from  the 
results  obtained  using  the  direct  simulation  Monte  Carlo  method. 

3  Some  authors  advocate  the  use  of  higher-order  continuum  equations  such  as  the  Burnett  equations 
(Ref.  6).  However,  there  is  an  inherent  problem  with  using  this  approach.  The  continuum  equations 
fail  because  the  velocity  distribution  function  becomes  highly  non-Maxwellian,  or  even  a  bi-modal 
mixture  of  perturbed  Maxwellians.  No  perturbation  to  a  single  Maxwellian  can  represent  such  a 
bi-modal  velocity  distribution  function.  Thus,  an  approach  based  on  higher-order  perturbations  to  a 
single  equilibrium  velocity  distribution  function  cannot  work. 
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2.2  Governing  Equations 

In  this  section  the  governing  equations  are  given;  further  information  on  their  deriva¬ 
tion  can  be  found  in  Refs.  13  and  29.  We  will  discuss  the  derivation  of  the  vibrational 
energy  equation  in  some  detail  because  it  provides  a  good  example  of  how  to  derive  con¬ 
servation  equations  for  nonequilibrium  flows.  Also  this  derivation  will  show  the  strong 
interaction  between  the  vibrational  state  and  the  chemical  reactions  in  a  nonequilibrium 
gas. 


2.2.1  Mass  Conservation 


The  mass  conservation  equation  for  species  s  is: 

dps  ,  d  (n  ^  ^  _ 

dt  dx  ■  +  PsVsi)  ~  Ws 5 


(2.2) 


where  again,  ps  is  the  species  density  and  ws  is  the  chemical  source  term.  The  mass- 
averaged  velocity  is  Uj  and  the  diffusion  velocity  of  species  s  is  vSj  in  the  Xj  direction.  The 
mass-averaged  velocity  is  obtained  using: 


Uj  — 


£ 


-u 


SJ  ) 


ns 

p= £p», 


(2.3) 

S  =  1  '  S=1 

where  ns  is  the  number  of  chemical  species.  The  diffusion  velocity  is  the  velocity  of  species 
s,  usj,  relative  to  the  mass-averaged  velocity: 


Vsj  —  'U'sj  U'j  • 


(2.4) 


2.2.2  Momentum  Conservation 


The  momentum  equation  has  the  familiar  form  except  for  the  presence  of  the  electric 


held,  Ep 


~  Ta)  =  YeZsNsEj, 

S=  1 


(2.5) 


where  eZs  is  the  charge  of  species  s  and  Ns  is  the  species  number  density.  The  pressure  p 
is  the  sum  of  the  partial  pressures: 


ns  ns  j~y 

=  (2-6) 
s=l  s=l  s 

where  R  is  the  universal  gas  constant  and  Ms  is  the  molecular  weight  of  species  s.  T  is 
the  translational  temperature  of  the  gas  mixture.  The  expression  for  the  shear  stress  will 
be  given  below. 
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2.2.3  Total  Energy  Conservation 

The  total  energy  conservation  equation  lias  the  form: 

dE  d  nS'  n,S' 

~dt  +  dx:^E  +  P")Uj  ~  TijUi  +  q°  +  ^jPsvsihs)  =  ^ ^eNsZsEiUi ,  (2.7) 

S=1  S=1 

where  E  is  the  total  energy  per  unit  volume,  q-j  is  the  total  heat  flux  vector,  and  hs  is  the 
species  s  specific  enthalpy.  These  quantities  will  be  discussed  in  more  detail  below. 


2.2.4  Vibrational  Energy  Conservation 

The  derivation  of  the  vibrational  energy  conservation  equation  is  non-trivial  because 
the  vibrational  state  is  coupled  to  the  chemical  state.  As  we  discussed  in  Section  1.2.2, 
molecules  that  are  highly  vibrationally  excited  are  more  likely  to  dissociate  than  the  average 
molecule.  Thus,  when  dissociation  occurs,  the  process  removes  more  than  the  average 
vibrational  energy  from  the  vibrational  energy  pool.  Likewise,  when  recombination  occurs 
the  newly  formed  molecule  may  be  formed  at  an  elevated  vibrational  level.  A  complete 
derivation  of  this  equation  is  available  in  Ref.  41. 

We  assume  that  there  is  a  single  diatomic  species  in  a  multi-species  gas  mixture.  In  a 
manner  similar  to  Clarke  and  McChesney  (Ref.  8),  we  let  fSa(xi,vSa,t)dxdvSa  represent- 
the  number  of  particles  in  vibrational  level  a  of  species  s  in  a  volume  dx  =  dx  {  dx 2  dx 3 
and  dvs  =  dvs  dvo  dvs  .  The  velocity  distribution  function  fs  is  defined  such  that 

0 a.  C,Q'^  ^ ex 2  0  ck  3  v  «/  oc 

/-Too 

fsa  dvSa  —  nSa,  (2.8) 

-OO 


where  nSa  is  the  number  density  of  molecules  in  level  a  of  species  s. 

Under  the  assumption  that  the  translational  energy  is  a  classical  energy  mode  and 
that  all  velocities  in  the  range  —  00  to  +00  are  allowed,  we  can  show  that  fSa  obeys  the 
Boltzmann  equation: 


dfs 


dt 


+  vs 


dfs 


dxi 


+  FS 


dfs 


H  dx.. 


=  ct  -  c: 


(2.9) 


where  FSa  is  the  external  force  per  unit  mass  acting  in  the  i  direction  on  the  particles, 
and  Cf  and  C~  represent  the  number  of  particles  created  and  destroyed  due  to  collisions 
of  species  s  in  level  a  per  unit  time  and  per  unit  phase  space  volume.  These  terms,  Cf 
and  C~  ,  are  the  collision  integrals. 
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The  conservation  of  vibrational  energy  equation  for  species  s  is  found  by  taking  the 
moment  of  the  Boltzmann  equation  with  respect  to  eSa  and  summing  over  all  vibrational 
energy  levels  a.  This  yields 


dfSa  , 

vSa.  —  dvSa 

1  dxi 

(Ct-C-)dv,„. 


(2.10) 


where  eSa  is  the  amount  of  vibrational  energy  per  molecule  of  species  s  in  level  a. 
The  convection  term  on  the  left  hand  side  of  (2.10)  can  be  manipulated  to  yield 


dqz 


d.  , 

“I”  Tj  (  Evs\V si  +  Ui) 


dx;  dx; 


(2.11) 


where  Evs  —  Ylansotesol  is  the  vibrational  energy  per  unit  volume,  and  qvst  = 
V  ns  es  us  is  the  vibrational  heat  flux  in  the  xs  direction.  If  the  external  forces  iT 
are  independent  of  the  velocity,  then  the  force  term  in  (2.10)  is  identically  zero  as  it  as¬ 
sumed  that  the  distribution  function  fSa  vanishes  at  infinity.  The  source  term  on  the  right 
hand  side  of  (2.10)  represents  the  rate  of  change  of  the  number  of  molecules  in  level  a  in 
species  s.  This  can  be  written  as 


/OO 

£»„  (c+-c- 

-oo 


dns 

dt  )  coil 


(2.12) 


Combining  the  above  expressions  we  have 
dK,„  d 


dt 


T  (-E'us  T  Evs  V si  T  Qvsi )  ^  ^  6Sa  ^  ^  ^  Qvib- 


(2.13) 


To  find  an  expression  for  QVib  we  formally  need  to  make  the  following  assumptions: 

1.  The  system  of  interest  is  a  dilute  mixture  of  vibrating-dissociating  molecules  and  atoms 
weakly  interacting  with  an  infinite  heat  bath. 

2.  The  Born-Oppcnheimer  approximation  holds  so  that  the  vibrational  states  are  uncou¬ 
pled  from  the  rotational  and  electronic  states  of  the  molecule. 

3.  The  interaction  Hamiltonian  which  causes  transition  between  vibrational  levels  can 
be  treated  as  a  perturbation  on  the  energy  of  the  vibrating  molecules.  Thus,  quantum 
mechanical  perturbation  theory  can  be  used  to  derive  the  master  relaxation  equation. 
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Under  these  assumptions  Heims  (Ref.  21)  showed 

=  ”aa  (Ns  -  2  ns)2  -  nSa  nSa  +  ^  ( aSianSi  -  aSainSa),  (2.14) 

i 

where  Ns  is  the  number  of  s  atoms,  ns  is  the  number  molecules  of  species  s,  aSlk  is  the 
transition  probability  per  unit  time  from  vibrational  level  l  to  k,  vSa  is  the  recombination 
rate  of  atoms  to  molecules  in  level  a,  and  fiSa  is  the  dissociation  rate  from  level  a. 

Evaluating  the  source  term  QVib  yields 


Qvib  vSa  (Ns  2  ns)  /iSa  nSa  +  ^  ^  (c 


n. 


n. 


(2.15) 


The  first  two  terms  of  (2.15)  are  respectively  the  gain  and  loss  of  vibrational  energy  due 
to  chemical  reactions.  The  last  two  terms  of  (2.15)  account  for  the  exchange  of  vibrational 
and  translational  energy  due  to  collisions.  At  the  microscopic  level,  these  two  processes 
are  not  linked  and  can  be  treated  independently.  Thus  we  define 

Qchem  —  ^  ^  £sa  (Ns  —  2  ns)  —  fiSa  ns (2-16) 

Oi 

and 

Qv-t  =  (°s-  “  a»«i  n*“)  •  (2-17) 

ex.  i 


To  find  an  expression  for  QChem  we  sum  (2.14)  over  all  a  levels  and  define 


7s 


E 


V* 


IN  =  Y, 

a 


k'Sg  IT-Sq 

n. s 


(2.18) 


where  is  the  recombination  rate  and  /is  is  the  dissociation  rate  of  s  molecules, 
gives 


This 


(2.19) 


where  ys  (Ns  —  2 ns)2  and  ns  fis  represent  the  rate  of  change  of  the  number  of  molecules 
of  species  s  per  unit  time  and  unit  volume  due  to  the  forward  and  backward  chemical 
reactions,  respectively.  In  (2.2),  let  ws  —  Wfs  +  Wbs  so  that  Wfa  and  Wf,a  are  the  rates  of 
change  of  the  mass  of  species  s  per  unit  time  and  unit  volume  due  to  the  forward  and 
backward  chemical  reactions,  respectively.  This  implies 


7S  (Ns  -  2 ns)2  =  —  ns  /is 
ma 


Wfs 

ms  ' 


(2.20) 
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The  average  vibrational  energy  gained  during  recombination  and  the  average  vibra¬ 
tional  energy  lost  during  dissociation  are  the  weighted  averages  of  the  level-specific  disso¬ 
ciation  rates  multiplied  by  the  vibrational  energy  of  each  level.  They  are  respectively: 


G 

E 


5Ea  7  s 

'Yhct  _  E^a  ^sa^sa^s 

SaAboTsa  ns 


(2.21) 


In  general,  the  nSa  in  the  above  equations  are  functions  of  T  and  Tv.  However  at 
equilibrium  T  —  Tv,  and  we  can  define  n*Q  as  the  number  density  at  equilibrium.  Thus 
we  have 

*»„  (V,  -  2n‘,f  -  ^  <)  =  0.  (2.22) 

i 

Now  we  make  the  assumption  that  at  equilibrium  each  process  must  be  in  equilibrium 
independent  of  the  other  process.  This  is  consistent  with  the  fact  that  we  are  treating 
Qchem  and  Q v _T  as  independent  on  a  microscopic  level.  This  gives 


ySa  (Ns  —  2  n*)2  —  fj,Sa  n*a  —  0.  (2.23) 

If  the  transition  probabilities  are  independent  of  time,  (2.22)  must  hold  for  all  time.  If 
we  substitute  the  expressions  for  v8a  from  the  above  equation  into  (2.16),  after  some 
manipulation  we  find 

EaVsan*Sa 

If  we  write  E  =  E(T,  Tv)1  we  immediately  see  that  G  =  E(T,T). 

Therefore,  any  physically  consistent  model  for  the  vibrational  energy  source  term  due 
to  chemical  reactions  must  be  of  the  form 


Qc^m  =  ——  \E(T,  Tv)  wfs  +E(T,T)wbs 

ms  \ 


(2.25) 


An  expression  for  QV_T  was  originally  derived  by  Landau  and  Teller  (Ref.  28)  for 
simple  harmonic  oscillators  not  undergoing  dissociation.  They  found 


^  E„.(T)  -  EV.(TV) 

Vv-T  — 

Tvib 

where  Tvtb  is  the  vibrational  relaxation  time,  and  is  given  by  a  theoretically  determined 
expression  as  a  function  of  the  local  thermodynamic  state  of  the  gas.  Under  more  general 


(2.26) 
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conditions,  QV_T  will  have  the  above  form,  but  tvh,  will  be  different  and  will  also  depend 
on  the  oscillator  model  used. 


The  equation  for  the  conservation  of  vibrational  energy  of  species  s  has  the  final  form: 


0EVS 

dt 


d 

“h  ^  —  (^Evs  Ui  +  Evs  V si  "l-  Qvsi )  Qv  —  t  T  Q 


(2.27) 


For  the  case  where  there  is  more  than  one  vibrationally  excited  species,  an  additional  term 
must  be  included  to  account  for  the  rate  of  vibrational  energy  transfer  to  species  s  from 
the  other  vibrationally  excited  species  (Ref.  50). 

If  the  vibrational  energy  modes  are  tightly  coupled,  then  there  will  be  a  single  vibra¬ 
tional  temperature  Tv,  and  the  total  vibrational  energy  equation  is: 

Q  _  _  _  _ 

~dt~  Ox  ^  l  EyaVsi  T  ^  ^  Qvsi )  y  ^  Qv-TS  T  y  ^  Q  ChemS-  (2.28) 

*  S  S  S  S 


where  Ev  =  Y,a  Evs- 


2.2.5  Additional  Internal  Energy  Conservation  Equations 

In  a  similar  fashion,  conservation  expressions  for  the  other  internal  energy  modes 
(rotational  and  electronic)  may  be  derived.  In  practice,  unless  the  conditions  of  interest  are 
at  very  low  density,  the  rotational  and  translational  energies  are  usually  considered  to  be  in 
equilibrium.  This  obviates  the  need  for  a  separate  rotational  conservation  equation.  Also, 
if  there  is  appreciable  population  of  the  excited  electronic  states,  it  is  often  assumed  that 
the  temperatures  that  characterize  the  free  electron  translational  energy  and  the  bound 
excited  electronic  energy  are  the  same.  That  leads  to  a  single  electron-electronic  energy 
conservation  equation.  In  many  cases,  it  is  valid  to  assume  that  the  vibrational  modes 
are  also  in  equilibrium  with  the  electron-electronic  energy  (because  of  resonant  coupling 
between  the  N2  vibration  and  free  electrons  (Refs.  30,  50)).  Then,  the  vibrational  energy 
equation  (2.13)  is  extended  to  include  these  effects. 

There  is  uncertainty  about  how  to  model  the  energy  in  the  free  electron  and  electronic 
energy  modes.  We  favor  the  approach  of  Gnoffo  et  al.  (Ref.  13)  in  which  it  is  assumed 
that  Tv  —  Te  —  Te£.  Namely,  that  because  of  the  strong  electron-vibration  coupling  in 
high  temperature  air,  we  typically  have  Tv  —  Te.  It  is  probably  reasonable  to  assume  that 
Te  =  Tej.  (There  really  is  no  other  alternative,  in  any  case.)  The  other  approach  is  to 
assume  that  T  —  Te  —  T,,( . 
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2.3  Electric  Field 

An  expression  for  the  electric  field  Et  may  be  derived  using  the  electron  momentum 
conservation  equation: 

d  d  ~ 

(Pe^h)T  Qrj.  T  Pe^ij')  &NeEi  A  Pei i  (2.29) 

where  we  have  neglected  the  mass  diffusion  and  shear  stress  terms.  Pet  represents  the 
momentum  transfer  between  the  electrons  and  the  heavy  particles  due  to  collisions.  For 
weakly  ionized  flows,  this  term  is  small  relative  to  the  electric  field  term.  Then,  if  we  take 
the  ratio  of  the  dynamic  pressure  of  the  electron  gas  to  the  electron  pressure  and  assume 
that  the  electron  speed  and  temperature  are  about  the  same  as  the  bulk  gas,  we  have 


PeVl^MeVl^M^^  Me 

pe  RTe  RT  M 


(2.30) 


Where  Ai  is  the  Mach  number,  and  M  is  the  average  molecular  weight  of  the  mixture. 
The  ratio  Me/M  is  of  the  order  of  10-6,  and  for  conditions  of  interest  the  square  of  the 
Mach  number  will  be  of  the  order  of  103  at  most.  Therefore,  we  can  neglect  the  electron 
dynamic  pressure  relative  to  the  electron  pressure,  and  the  steady-state  electric  field  may 
be  expressed  as 


1  dpe 
Nee  dxi 


(2.31) 


This  expression  for  the  electric  field  may  be  inserted  in  the  momentum  equation,  (2.5), 
and  the  total  energy  equation,  (2.7).  Generally,  this  term  has  little  effect  on  the  flow  field 
of  atmospheric  entry  vehicles  and  is  often  neglected. 
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2.4  Equations  of  State 

The  relationship  between  the  conserved  quantities  and  the  non-conserved  quantities 
such  as  pressure  and  temperature  are  discussed  in  this  section.  The  total  energy,  E,  is 
made  up  of  the  separate  components  of  energy: 

ns  ns  ns  ns  n 

E  =  Ets  +  Ers  +  Evs  +  Ee  +  Eeis  +  \puiUi  +  Psh°s ,  (2.32) 

S^e  s#e  s^e  s#e  s^e 

which  are  the  translational,  rotational,  vibrational,  electron  translational,  electronic,  ki¬ 
netic,  and  chemical  energies,  respectively.  The  heavy  particle  translational  energy  is 
Ets  —  PsCvtsT ,  where  cvst  —  § R/Ms  and  T  is  the  translational  temperature.  The  ro¬ 
tational  energy  is  Ers  —  pscvrsTri  and  cvrs  —  R/Ms  if  the  particle  has  two  degrees  of 
rotational  freedom;  Tr  is  the  rotational  temperature.  As  discussed  above,  for  the  case 
where  the  rotational  energy  modes  relax  quickly,  we  can  assume  that  the  rotational  en¬ 
ergy  modes  are  equilibrated  with  the  translational  energy  modes.  Then,  we  can  write 
Ets  +  Ers  —  PsCvsT ,  where  cvs  and  T  are  the  translational-rotational  specific  heat  and 
temperature,  respectively. 


The  vibrational  temperature  of  species  s  is  determined  by  inverting  the  expression  for 
the  vibrational  energy  contained  in  a  simple  harmonic  oscillator  at  the  temperature  Tvs: 


E 


Vs 


ps^vs  -  Ps 


R  Oys 

]\/[g 


-1’ 


(2.33) 


where  9VS  is  the  characteristic  temperature  of  vibration.  If  we  assume  that  there  is  a  single 
vibrational  temperature,  we  must  invert  a  more  complicated  expression  for  Tv 


Ey  —  ^  ^  Ps 
s 


R  9ys 

Ms  edvs/Tv  —  1 


(2.34) 


The  free  electron  translational  energy  is  given  by  Ee  —  pecveTe,  where  cve  =  § R/Me. 
In  cases  where  the  electron  energy  is  assumed  to  be  in  equilibrium  with  another  energy 
mode,  a  different  temperature  is  used  in  the  above  expression. 

As  discussed  above,  the  total  pressure  is  the  sum  of  the  partial  pressures, 

ns  „ 

p  =  'l2psM~T +Pe’  (2-35) 

s^e 

and  the  electron  pressure  is  given  by 

Pe  =  peWTe'  ('2'36'1 
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The  enthalpy  per  unit  mass,  hs,  is  defined  to  be 

ha  —  cvtsT  +  cvrsTr  +  evs  +  ee£s  +  h°  +  — .  (2.37) 

Ps 

The  expression  for  the  energy  contained  in  the  excited  electron  states  comes  from  the 
assumption  that  they  are  populated  according  to  a  Boltzmann  distribution  governed  by 
the  electronic  temperature,  Te£.  This  yields 

_  R  Pis^eiis  exP(~^elas/^el)  [r,  oo\ 

&els  i  i  v— \oo  /  a  /rri  )  (y.oo) 

Ms  Li=0  9ia  expi-Vetis/Tet) 

where  giS  is  the  degeneracy  of  the  excited  electronic  state  i  and  9e£is  is  the  excitation 
energy  of  that  state  (Ref.  50).  Usually  only  the  first  several  terms  in  these  summations 
are  required  for  hypersonic  applications. 

Another  approach  is  to  use  the  fits  for  pure  species  thermodynamics  data  from  Gordon 
and  McBride  (Ref.  15).  With  this  approach,  the  species  translational- rotational  energy  and 
chemical  energy  can  be  subtracted  from  the  total  species  energy  to  obtain  the  vibrational- 
electronic  energy.  These  curve-fits  are  valid  to  high  temperature  (20,  000  K)  and  result  in 
a  more  accurate  characterization  of  the  internal  energy. 


2.5  Diffusion  Velocity,  Shear  Stress,  and  Heat  Flux 

The  shear  stresses  are  assumed  to  be  proportional  to  the  first  derivative  of  the  mass- 
averaged  velocities,  and  the  Stokes  assumption  for  the  bulk  viscosity  is  made.4 

Therefore  the  expression  for  the  shear  stress  tensor  is 


Tij  ~  Py 


( dui  du 


+ 


dxj  dx 


±\dukz 
+  A  — —  o 

ox  I 


IJ , 


x=~r- 


(2.39) 


And  the  heat  conduction  vectors  are  assumed  to  be  given  by  the  Fourier  heat  law 


dT 


Qtrj  —  ^ 


dxj: 


dT 

u ^  vs 

Qvsj  —  ^vs~7\  • 

dxj 


(2.40) 


4  The  subject  of  bulk  viscosity  is  interesting.  It  has  long  been  recognized  that  the  internal  energy 
modes  of  a  gas  may  affect  the  speed  of  sound  (Refs.  18,  28,  57).  This  is  commonly  attributed  to  the 
bulk  viscosity  because  sound  propagation  is  a  dilatational  process,  and  the  bulk  viscosity  provides 
a  means  of  changing  the  speed  of  sound  due  to  non-zero  dilatation.  However,  this  is  not  physically 
consistent.  Instead,  the  proper  finite-rate  internal  energy  relaxation  equations  should  be  used  to 
obtain  the  correct  speed  of  sound. 
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In  some  cases,  it  is  preferable  to  represent  the  vibrational  energy  conduction  with  the 
gradient  of  the  specific  vibrational  energy. 

Qf> 

Qvsj  ^vs~^  • 

OXj 

Except  at  very  high  enthalpy,  we  find  that  the  vibrational  energy  gradients  give  more 
robust  results. 

A  viscosity  model  for  reacting  air  developed  by  Blottner  et  al  (Ref.  3)  may  be  used  to 
determine  the  species  viscosity,  fj,s.  This  work  uses  kinetic  theory  to  find  curve-fit  expres¬ 
sions  for  the  viscosity  of  each  species.  Another  source  of  expressions  for  nonequilibrium 
calculations  may  be  found  in  the  work  of  Gupta  et  al.  (Ref.  19).  This  information  is  too 
extensive  to  include  in  these  notes.  Excellent  reviews  of  the  calculation  of  species  transport 
properties  are  given  by  Palmer  and  Wright  (Refs.  43,  44). 

The  approximate  conductivity  of  the  translational-rotational  and  vibrational  temper¬ 
atures  for  each  species  may  be  derived  from  an  Eucken  relation  (Ref.  62).  With  this 
approach,  it  is  assumed  that  the  transport  of  translational  energy  involves  correlation 
with  the  velocity,  but  the  transport  of  internal  energy  (rotational  and  vibrational)  has  no 
correlation.  The  result  is  that 


KS 


l^s  (  2  Cuts  T  Ct,7-S)  5 


^vs  Vv  l^s^vvsj 


(2.42) 


where  r)v  —  1.2  is  derived  from  kinetic  theory  (Ref.  42),  and  cvvs  is  the  species  s  vibra¬ 
tional  specific  heat.  When  vibrational  energy  derivatives  are  used,  the  transport  coefficient 
becomes: 

kvs  —  T]vij,s.  (2.43) 

This  approach  based  on  the  Eucken  relation  is  approximate  and  is  only  valid  up  to  about 
6000  K  for  air;  at  higher  temperatures  a  more  sophisticated  approach  must  be  used.  Palmer 
and  Wright  (Refs.  43,  44)  provide  a  quantitative  assessment  of  the  available  approaches. 

Once  the  pure  species  viscosity  and  conductivity  have  been  computed,  the  mixture 
properties  must  be  obtained.  This  is  often  done  with  the  Wilke  semi-empirical  mixing  rule 
(Ref.  63),  however  Palmer  and  Wright  show  that  this  approach  is  subject  to  serious  error. 
They  recommend  the  use  of  the  Armaly-Sutton  (Ref.  1)  mixing  rule  because  it  is  more 
accurate  and  less  costly  than  the  solution  of  the  full  multi-component  diffusion  equations. 
However,  the  parameters  in  the  Armaly-Sutton  model  may  need  to  be  tuned  for  particular 
gas  mixtures  and  conditions  (Refs.  43,  44). 
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If  we  assume  that  the  diffusive  fluxes  due  to  pressure  and  temperature  gradients  are 
negligible,  then  the  diffusion  velocity  of  each  component  of  the  gas  mixture  is  proportional 
to  the  gradient  of  the  mass  fraction.  With  the  additional  assumption  of  binary  diffusion 
where  species  s  diffuses  into  a  mixture  of  similar  particles,  we  have 

Qc 

psvsj  =  -pDs-^~.  (2.44) 


The  diffusion  coefficient,  DS:  is  derived  by  assuming  a  constant  Lewis  number,  Le,  which 
by  definition  is  given  by 


Le  — 


pDCf 

K 


(2.45) 


For  air,  Le  is  typically  taken  to  have  a  value  of  1.4,  and  thus  the  uncharged  particles 
all  have  the  same  D,  but  the  diffusion  coefficient  for  ions  is  assumed  to  be  doubled  (the 
ambipolar  diffusion  assumption  holds)  because  of  the  existence  of  an  electric  field. 


A  much  more  accurate  approach  for  computing  diffusion  in  a  gas  mixture  is  the  self- 
consistent  effective  binary  diffusion  (SCEBD)  approach  of  Ramshaw  and  Chang  (Refs.  54, 
55)  and  the  extension  of  the  method  to  multi-temperature  plasmas  (Ref.  56).  This  ap¬ 
proach  has  been  shown  to  yield  accurate  results  for  high-enthalpy  atmospheric  entry  flows 
(Ref.  16).  At  atmospheric  entry  conditions  (particularly  when  the  gas  is  ionized),  the 
constant  Lewis  number  approximation  is  not  valid,  and  the  SCEBD  approach  should  be 
used. 


2.6  Internal  Energy  Relaxation  Rates 

The  rate  of  energy  exchange  between  vibrational  and  translational  modes  has  been 
discussed  extensively  (Ref.  28).  The  rate  of  change  in  the  population  of  the  vibrational 
states  at  low  temperatures  is  described  well  by  the  Landau- Teller  formulation  where  it  is 
assumed  that  the  vibrational  level  of  a  molecule  can  change  by  only  one  quantum  level  at 
a  time.  The  resulting  energy  exchange  rate  is 


Qv-i 


AT) 


. '■■■  Ts  L_T  '  V 


(2.46) 


Where  e*a{T )  is  the  vibrational  energy  per  unit  mass  of  species  s  evaluated  at  the  local 
translational  temperature  and  <  tSl_t  >  is  the  molar  averaged  Landau- Teller  relaxation 
time 


<r. 


SL-T 


>= 


ErXr/T, 


for  r^e. 


sr  i^_T 


(2.47) 
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An  expression  developed  by  Millikan  and  White  (Ref.  39)  yields  the  Landau- Teller  inter¬ 
species  relaxation  times,  rsrL_T,  in  seconds  using  the  function 

t sri^_T  —  -  expL4sr(T-1/3  —  0.015/ig/4)  —  18.42],  p  in  atm, 

p 

Asr  =  1.16  x  lO“3/x^204/3,  (2‘48) 

Psr  =  MsMr/(Ms  +  Mr). 


There  are  notable  exceptions  to  the  Millikan- White  formula,  particularly  for  N2  and 
O2  vibration-translation  relaxation  involving  atomic  oxygen,  and  the  relaxation  of  CO2 
(Ref.  50). 

A  modification  to  the  translational-vibrational  relaxation  rate  is  made  to  account  for 
the  limiting  collision  cross-section  at  high  temperatures.  The  Landau- Teller  rate  expres¬ 
sion  from  Millikan  and  White  yields  a  relaxation  rate  that  is  unrealistically  large  at  high 
temperatures  due  to  an  overprediction  of  the  collision  cross-section.  The  addition  of  the 
limiting  cross-section  rate  corrects  this  inaccuracy.  As  suggested  by  Park  (Ref.  50),  a  new 
relaxation  time,  rvs,  that  is  the  sum  of  the  Landau- Teller  relaxation  time  and  the  collision- 
limited  relaxation  time,  rcs,  corrects  this  inadequacy.  Thus  if  we  use  (2.46)  with  this  new 
rate,  we  have  the  final  form  of  the  translational-vibrational  energy  exchange  rate: 


Q 


V  —  T  s 

Tvs 


e*vs(T)  -  ^ 


<  Tsl_t  ^  A  Tcsi 


(2.47) 


where 


7~r.s  — 


c.a. 


VNS 


(2.49) 


cs  is  the  average  molecular  speed  of  species  s,  cs  =  \J&RT /tiMs,  and  Ns  is  the  number 
density  of  the  colliding  particles.  The  expression  for  the  limiting  collision  cross-section, 
av .  is  assumed  to  be  as  given  by  Ref.  50: 


av  =  1(T17(50,000/T)2  cm2, 


(2,51) 


where  T  is  in  K.  This  expression  was  originally  developed  for  nitrogen,  but  has  been  applied 
to  the  other  diatomic  molecules. 


2.7  Chemical  Source  Terms 

The  source  term  for  each  chemical  species  may  be  constructed  using  the  law  of  mass 
action  (Ref.  62)  and  a  given  set  of  chemical  reactions.  In  this  section,  we  develop  the 
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chemical  source  terms  for  a  simple  chemical  kinetics  model;  it  is  easy  to  generalize  this 
discussion  to  other  models. 

For  high  temperature  non-ionized  air  there  are  five  primary  components,  which  may 
be  ordered  as  follows,  N2,  O2,  NO,  N,  and  O.  The  most  important  chemical  reactions 
between  these  species  are 

N2  +  M  ^  2N  +  M 
02  +  M^20  +  M 

NO  +  M^N  +  O  +  M  (2.52) 

N2  +  O  ^  NO  +  N 
NO  +  O  ^  02  +  N, 

where  M  represents  any  particle  that  acts  as  a  collision  partner  in  the  reaction.  The 
first  three  are  dissociation  reactions  and  the  remaining  two  are  exchange  reactions.  Each 
reaction  is  governed  by  forward  and  backward  reaction  rate  coefficients,  kfm  and  kbrn , 
respectively.  These  five  reactions  may  be  written  in  order  in  terms  of  the  reaction  rates  as 


fti  = 


= 


n 

m 

£[- 


it  1 
J  1  r 


'fir 


+  hlr 


Pn  Pn  P n 
Mn  Mn  M, 


P  N2  Pm 
Mn2  Mm 

P°2  Pm  I  kh _ 

Mq2  Mm  b2m  M0  M0  Mri 


P O  P O  Pm 


^3  = 
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72-4  =  —  kf4 
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Mn2  Mo 
Pno  Po 
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+ 


Mno  Mn 
Po2  Pn 
Mq  2  Mn 


(2.53) 


Thus,  the  source  terms  that  represent  the  inter-species  mass  transfer  rates  may  be  con¬ 
structed  as 

^n2  —  Mn2  (72-i  +  72.4) 
wo2  —  Mo2  (n2  —  72.5) 

wno  —  7V^no(7^3  —  72-4  +  72-s)  (2-54) 

wn  =  Mn( — 272-1  ~~  72-3  —  72-4  —  72-5 ) 

wo  —  Mo(— 272. 2  —  72-3  +  72-4  +  72-5 ) . 

We  should  note  that  the  sum  of  the  mass  transfer  rates  is  identically  zero  and  that  elemental 
conservation  holds,  as  required. 
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In  equilibrium,  the  forward  and  backward  reaction  rates  of  reaction  rn  have  the  func¬ 
tional  form: 


kfj?) 
h  JT) 


CfrnTVrn  exp(— 0m/T), 

kfJT ) 

^JT)  ’ 


(2.55) 


where  the  constants  Cfm,  rjm ,  and  6m  are  experimentally  determined  (e.g.  Refs.  51,  52) 
or  computed  using  computational  chemistry,  and  Keqm  is  computed  from  first  principles 
using  thermodynamic  data  (Ref.  15).  However  as  discussed  by  many  authors  and  as  shown 
in  Sec.  2.2.4,  the  vibrational  state  of  the  gas  affects  the  dissociation  rate.  Many  models  for 
the  vibration-dissociation  coupling  process  have  been  proposed.  We  will  discuss  the  most 
widely  used  of  these  models  here. 


The  most  widely  used,  because  of  its  simplicity,  is  the  Park  TTV  model,  in  which  the 
temperature  that  governs  the  forward  reaction  rate  is  replaced  by  an  effective  or  average 
temperature,  Ta.  Park  originally  proposed  that  Ta  —  \jTTy ,  however  a  more  appropriate 
expression  is 

Ta  =  T^T^,  (2.56) 

where  cf)  is  usually  taken  as  0.7.  This  model  is  based  on  some  more-or-less  heuristic  reason¬ 
ing,  but  it  seems  to  work  well  and  gives  reasonable  results.  However  as  we  saw  during  the 
derivation  of  the  vibrational  energy  conservation  equation,  the  chemical  reaction  rate  im¬ 
plies  a  certain  rate  of  vibrational  energy  removal  due  to  the  reactions.  This  energy  removal 
rate  is  state-specific,  so  unless  the  reaction  rate  model  is  state  specific,  it  is  impossible  to 
derive  an  appropriate  QChem  (see  (2.25)).  Typically,  it  is  assumed  that  QChem  —  0.3 De, 
where  De  is  the  dissociation  energy. 


Other  authors  have  used  a  more  detailed  derivation  of  the  vibration-dissociation  cou¬ 
pling  process.  For  example,  the  CVDV  (coupled  vibration-dissociation-vibration)  model  of 
Marrone  and  Treanor  (Refs.  37,  61)  assumes  a  Boltzmann  distribution  of  the  vibrational 
states,  and  allows  preferential  removal  due  to  dissociation  from  the  upper  states.  This  re¬ 
sults  in  an  effective  dissociation  rate  that  is  a  function  of  the  vibrational  and  translational 
temperatures  and  the  parameter  U: 


kf  = 


Q(T)Q(Tr)  „/T 

Q(TV)Q(-U)  > 


(2.57) 


where  Q(T)  is  the  vibrational  partition  function  evaluated  at  temperature  T  and  is  defined 
as 

N 


Q(T)  =  J2e~6a/kT- 

a=0 


(2.58) 
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In  (2.57)  TF  represents  a  modified  temperature  and  is  given  by 

1  _  i  1  1 

T~F~  Yv~  T~  U' 


(2.59) 


For  U  —  oo  there  is  an  equal  probability  of  dissociation  from  all  levels,  and  as  U  decreases, 
the  probability  increases  that  a  dissociating  molecule  comes  from  an  upper  vibrational 
energy  level.  U  —  1/3  is  typically  used. 

The  CVDV  model  naturally  results  in  an  expression  for  QChem ■  The  expressions  for 
the  energy  removed  and  gained  during  reactions  found  in  (2.25)  are: 


E(T,TV) 
E(T ,  T) 


1  N 
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€(y  C 


—ea/kTF 


^~U)  to 
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en  e 


Y/kU 


(2.60) 


Knab,  Fruhauf  et  al.  (Refs.  25-27)  developed  the  CVCV  (coupled  vibration-chemistry- 
vibration)  model  that  generalizes  the  approach  of  Marrone  and  Treanor.  This  leads  to 
expressions  for  the  effective  reaction  rate  and  the  vibrational  energy  loss  terms  that  are 
similar  to  the  CVDV  results. 

The  Macheret  and  Rich  (Ref.  36)  model  takes  a  classical  approach  to  the  coupling 
problem.  As  opposed  to  the  discrete  energy  levels  of  the  real  oscillator,  Macheret  and  Rich 
assume  that  the  vibrational  energy  distribution  function  can  be  approximated  by 


/(O 


kTv 
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kTv 
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ev/kTv 


e  1  _  tu-q 
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if  <  £i 


if  ev  >  ei, 


(2.61) 


where  e  \  is  approximately  one  half  of  the  dissociation  energy. 

This  distribution  function  is  an  attempt  to  take  into  account  the  fact  that  the  vibra¬ 
tional  energy  mode  does  not  relax  through  a  series  of  Boltzmann  distributions.  Macheret 
and  Rich  assume  that  the  nonequilibrium  distribution  function  can  be  characterized  by 
a  Boltzmann  distribution  at  temperature  Tv  for  the  lower  levels  and  another  Boltzmann 
distribution  at  temperature  T  for  the  upper  levels.  The  expression  for  vibrational  energy 
becomes 

p  (rp  rp\  _  PN2  Jo  evf(tv)dev 

TON2  f0  f(ev)dev 
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where  De  is  the  dissociation  energy  of  the  molecule.  The  vibrational  energy  is  now  a  weak 
function  of  the  translational-rotational  temperature. 

Macheret  and  Rich  generalize  the  Arrhenius  formula  for  vibrational  nonequilibrium 
by  considering  a  threshold  energy  function  which  determines  the  minimum  total  energy 
in  a  collision  necessary  for  dissociation.  The  concept  of  preferential  removal  is  built  into 
this  method  by  the  theoretically  determined  threshold  function.  This  method  also  ac¬ 
counts  for  the  rotational  state  of  the  molecule  and  can  be  used  for  flows  with  rotational 
nonequilibrium. 

The  uonequilibrium  dissociation  rate  is  found  to  be 

kf  =  k0{ki  +  hi  +  kh),  (2.63) 

where  kg,  kt .  and  kh.  are  the  rates  from  the  low,  intermediate,  and  high  vibrational  levels. 
These  expressions  are  complicated  and  will  not  be  repeated  here  (Ref.  36).  This  approach 
also  yields  an  expression  for  QChf,„,  in  terms  of  the  parameters  of  the  model  and  the 
vibrational  and  translational  temperatures. 

There  is  one  critical  issue  associated  with  the  use  of  these  vibration-dissociation  mod¬ 
els.  It  is  sometimes  difficult  to  interpret  the  experimental  data  used  to  derive  the  constants 
in  the  Arrhenius  expression  for  the  forward  reaction  rate  (2.55).  In  many  cases  the  reac¬ 
tion  rates  were  measured  in  shock-heated  gas  when  the  gas  may  be  in  thermo-chemical 
nonequilibrium.  In  this  case,  it  is  important  to  interpret  the  experimental  data  in  a  man¬ 
ner  that  is  consistent  with  the  vibration-dissociation  model  being  used.  For  example,  Park 
(Refs.  48,  49)  made  an  extensive  study  of  air  reaction  rates  in  the  light  of  the  TTV  model. 

The  modeling  of  vibration-dissociation  coupling  is  still  an  open  issue,  and  virtually  no 
work  has  been  done  in  this  area  in  the  past  15  years.  It  may  be  that  with  recent  results 
from  computational  chemistry,  it  will  be  possible  to  study  the  dissociation  process  in  much 
greater  detail  and  many  of  these  issues  will  be  resolved. 

2.8  Boundary  Conditions 

The  boundary  conditions  for  hypersonic  flows  can  range  from  very  simple  (isothermal, 
non-catalytic  surface)  to  extremely  complicated  (mass  injection  with  in-depth  material 
response).  In  this  section,  we  cover  only  a  few  of  the  simpler  boundary  conditions. 

Usually  it  is  appropriate  to  assume  that  there  is  no  slip  at  the  body  surface,  and 
therefore  the  velocity  on  the  surface  is  zero.  Often,  the  wall  temperature  is  either  specified 
due  to  material  properties  or  the  mode  of  operation  of  a  particular  test  facility.  Seldom 
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is  an  adiabatic  wall  condition  used  because  at  hyp ervelo city  conditions  this  results  in 
unrealistically  high  surface  temperatures.  There  are  several  situations  that  require  more 
complicated  surface  boundary  conditions. 

At  low  densities,  there  may  be  velocity  and  temperature  slip  at  the  wall.  That  is,  if 
the  Knudsen  layer  at  the  wall  has  appreciable  thickness,  the  velocity  at  the  vehicle  surface 
may  not  approach  zero.  Gokcen  (Ref.  14)  and  others  have  developed  expressions  that  may 
be  used  to  calculate  the  velocity  and  temperature/energy  slip. 

When  the  surface  promotes  recombination  of  the  gas,  a  finite-rate  wall  catalysis  model 
must  be  used.  Typically,  the  wall  catalysis  is  expressed  as  a  catalytic  efficiency  of  a  surface 
reaction,  ar,  in  the  expression: 


(2.64) 


where  Ms  is  the  molecular  weight  of  the  species  that  is  recombining  at  the  wall,  and 
Tw  is  the  wall  temperature.  ar  is  measured  experimentally,  and  is  often  a  function  of 
temperature.  Then,  the  mass  flux  of  the  recombined  species  at  the  surface  is  pskr.  The 
relevant  boundary  condition  for  the  surface  state  can  be  obtained  by  equating  this  mass 
flux  to  the  diffusive  mass  flux  of  the  species  given  by: 


(2.65) 


Using  the  assumption  that  the  normal-direction  pressure  gradient  is  zero  and  the  boundary 
condition  for  temperature,  the  state  of  the  gas  on  the  wall  may  be  computed  using  an 
iteration  scheme. 

In  some  cases,  the  vehicle  may  fly  at  a  free-stream  condition  for  long  enough  that 
the  surface  reaches  a  locally-constant  temperature.  When  there  is  no  re-radiation  from 
the  surface,  this  is  the  adiabatic  wall  condition.  However,  in  cases  where  there  is  a  high 
surface  temperature,  surface  re-radiation  is  important.  Then,  the  convective  heat  transfer 
to  the  surface  is  balanced  by  the  black-body  re-radiation  heat  transfer  rate,  qrad  =  fJcT,),. 

At  very  high  heating  rates,  the  surface  may  ablate.  Then,  the  processes  of  oxidation, 
sublimation,  and  spallation  must  be  considered.  These  processes  are  complicated  and  are 
beyond  the  scope  of  these  notes. 
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3.  Numerical  Methods  for  Hypersonic  Flows 

In  this  section  we  discuss  numerical  methods  that  are  appropriate  for  solving  the 
governing  equations  discussed  above.  We  will  focus  on  one  method,  data-parallel  line- 
relaxation  (DPLR)  (Ref.  67)  with  modified  Steger- Warming  flux  vector  splitting  (Ref.  60). 
This  approach  has  direct  connection  to  other  more  modern  upwind  methods,  and  has  been 
shown  to  be  reliable  for  a  wide  range  of  applications.  As  we  described  above,  the  method 
used  must  be  parallelizable  and  implicit  so  that  solutions  may  be  obtained  in  a  reasonable 
amount  of  time. 

We  consider  a  gas  composed  of  ns  species  whose  translational  and  rotational  modes 
are  in  equilibrium.  We  assume  that  the  vibrational  state  is  characterized  by  a  single 
vibrational  temperature  because  of  strong  vibration-vibration  coupling.  Also,  there  is  no 
ionization.  This  gas  model  can  be  generalized  to  include  other  effects. 


3.1  Conservation-Law  Form  of  the  Governing  Equations 


The  governing  equations  for  the  nonequilibrium  flow  that  were  presented  in  the  pre¬ 
vious  section  may  be  written  in  a  form  that  is  more  suitable  for  the  derivation  of  the 
numerical  method.  This  is  the  conservation-law  form  of  the  differential  equations  where 
the  time  rate  of  change  of  the  vector  of  conserved  quantities  is  balanced  by  the  gradients  in 
the  flux  vectors  and  the  source  vector.  In  two  dimensions  the  governing  equations  written 
in  this  form  are 


du_  ar  og  _ 

dt  dx  dy 

where  the  vector  of  conserved  quantities,  U,  is  given  by 


(3.1) 


U  =  (pi,  p2, 


Prisi  pu ,  pv,  E,,.  E) 


(3.2) 


The  quantities  u  and  v  are  the  mass-averaged  velocity  components  in  the  x  and  y  directions 
respectively.  The  x  and  y  direction  fluxes  are  written  as 


/ 


pi(u  +  ui)  \ 

p2(u  +  u2) 


F  = 


\(E  +  p 


Pns  (w  T  Uns) 

PU2  +  p~  Txx 
plLV  7~Xy 

EyU  +  o  V/gEyg  +  Qvx 
'Txx'j'U'  7~xyV  +  qt  rx  +  q  vx 


(3.3) 
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Pl(v  +  Vi) 
P2(V  +  V2) 


\ 


G  = 


\(E  +  p 


Pnsiy  “I"  Vns) 
plLV  7~yX 
pV2  +P~Tyy 

EyV  +  VsEvs  -|-  Qvy 

7~yy)v  Tyx 'll  “h  Qtry  H“  Qvy  H-  Ps^s^s  ' 


(3.4) 


where  the  quantities  us  and  vs  are  the  x  and  y  components  of  the  diffusion  velocity  of 
species  s.  The  source  vector  is  made  up  of  terms  that  represent  the  mass,  momentum,  and 
energy  transfer  rates,  and  may  be  written  as: 


/  Wl  \ 

W2 


W 


Wns 

0 

0 


+  ESQ 
o 


chemS 

/ 


(3.5) 


3.2  An  Implicit  Finite- Volume  Method 

In  two  dimensions,  the  hnite-volume  approach  discretizes  the  flowheld  on  a  grid  of 
triangular  or  quadrilateral  elements.  The  x  and  y  locations  of  the  volume  corners  (nodes) 
are  stored,  and  the  state  of  the  gas  is  represented  with  volume-averaged  quantities  stored 
at  the  centroids  of  the  elements.  Each  face  has  an  outward-pointing  surface  normal  vector, 
S,  and  the  volume  of  each  element  is  given  by  V.  Figure  3.1  gives  a  graphical  representation 
of  this  scheme.  It  is  easy  to  extend  this  representation  to  three  dimensions. 
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Fig.  3.1  Numbering  scheme  for  the  two-dimensional  finite-volume  method. 


If  we  integrate  the  conservation  equations  over  a  finite  volume  cell,  V,  we  obtain: 

4-  [  UdV+  [  V  ■  F  dV  —  [  WdV  (3.6) 

ot  Jv  Jv  Jv 

dJL  +  lJsP.dS  =  W  (3.7) 

Where  F  —  FT+  Gj,  and  S  is  the  cell  surface- normal  vector.  For  a  finite  volume,  we  can 
then  interpret  this  expression  as 


dU 

~dt 


+ 


faces 


w. 


(3.8) 


where  we  have  dropped  the  bars,  and  we  sum  the  fluxes  across  the  faces  of  the  finite 
volume. 

Now,  we  can  obtain  a  fully  implicit  method  by  evaluating  the  fluxes  and  the  source 
vector  at  the  future  time  level,  n  +  1: 

un+i_un  +  ^_  J2Fn+1-S  =  Atwn+1.  (3.9) 

faces 

And  then  we  can  linearize  the  fluxes  and  source  vector  as: 

F)F 

Fn+1  pn  +  ^ / jjn+l  _  Jjn\  _  pn  _|_  j^n  / jjn+ 1  _  jjn\ 

(3-10) 

wn+1  ~  Wn  +  (un+1  -  un)  =  wn  +  cn  (un+1  -  un ) 
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Where  A  —  +  =  Az  +  Bf.  If  we  define  6Un  —  Un+1  —  Un,  we  obtain  the  expression: 


5Un  +  —  V  An5Ur 

V  LL 


■S-AtCn5Un  =  V  Fr 


S  + At  Wr 


(3.11) 


We  need  to  determine  how  to  evaluate  the  fluxes  at  the  cell  surfaces,  given  the  flow 
quantities  at  the  cell  centroids. 

3.2.1  Flux- Vector  Splitting 

Consider  only  the  inviscid  portion  of  the  fluxes;  the  viscous  fluxes  are  diffusive  and 
it  is  relatively  straight-forward  to  evaluate  them.  First,  let’s  derive  a  simple  first-order 
accurate  method  developed  by  Steger  and  Warming  (Ref.  60).  For  convenience,  define  a 
rotated  flux  vector  F'  such  that 

F'  —  Fsx  +  Gsy  (3.12) 

where  sx  and  sy  are  the  direction  cosines  of  the  surface  normal  vector,  S. 

We  recognize  that  the  inviscid  part  of  F'  is  homogeneous  in  U,  and  therefore  F'  =  A'U. 
We  can  diagonalize  A'  and  split  the  eigenvalues  into  those  that  are  positive  and  those  that 
are  negative: 

F  —  A'U  —  X~1AA'X  U 

=  X~xK\,X  U  +  X-1A^,X  U  (3.13) 

=  A’+U  +  A'~U  =  F,+  +  F'~ 

Where  is  the  diagonal  matrix  of  the  eigenvalues  that  are  positive,  and  contains 
the  negative  eigenvalues.  Physically,  F,+  represents  the  flux  moving  in  the  surface-normal 
direction,  and  F'~  is  the  flux  moving  in  the  opposite  direction.  Therefore,  when  we  evaluate 
the  fluxes  at  a  cell  face,  we  should  use  information  taken  from  the  appropriate  location,  as 
seen  in  Fig.  3.2. 


Fig.  3.2  -  Illustration  of  fluxes  across  a  surface. 
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Thus,  to  compute  F'  at  the  surface  shown  above,  we  use 


rUi.^Fu; 


(3.14) 


In  this  expression,  Steger  and  Warming  evaluated  the  Jacobians  at  the  same  location  as 
the  solution  vector;  however  this  leads  to  a  very  dissipative  method  (Ref.  34).  Instead,  it 
is  much  better  to  average  across  the  face  to  evaluate  the  Jacobians  using: 


ry..i  =  4i  ivr+A'-.u, 


i+ 


i+4 


i+ 1- 


(3.15) 


where 

A'±h=A,±(±(Ui  +  Ui+1)),  (3.16) 

that  is,  we  use  the  average  of  the  flow  quantities  on  either  side  of  the  cell  surface  to  evaluate 
the  Jacobian.  This  results  in  much  less  dissipation  than  the  original  Steger- Warming 
method. 

Interestingly,  the  above  approximation  for  the  flux  can  be  written  in  a  different  form 
by  combining  the  flux  components: 

F'+ 1  =  F\\(Ui  +  Ui+ 1))  -  l(X-1\\A,\X)(Ui+1  -  Ui)  (3.17) 

2  £ 


where  the  first  term  is  an  unbiased  average  of  the  flux  at  the  surface,  and  the  second  term 
is  an  upwind-biased  dissipative  flux.  This  is  the  familiar  Roe  form  of  the  flux. 


This  expression  can  be  used  to  obtain  a  first-order  accurate  approximation  to  the  flux 
at  each  face  of  the  element.  However,  for  strong  shock  waves  and  other  discontinuities, 
the  unbiased  averaging  across  the  face  will  produce  aphysical  results  (negative  densities 
and  energies).  Therefore,  a  sensor  must  be  used  to  smoothly  switch  back  to  the  more 
dissipative  form  of  the  flux  in  regions  of  strong  pressure  gradient.  We  use  a  weight  of  the 
form: 


w  —  1  — 


1  1 

2Kf  +  i! 


Pi  ~  Pi+i 
min(pj,pi+i) 


(3.18) 


Where  w  and  1  —  w  are  used  to  weight  Ui  and  Ui+ i  so  that  the  flux  is  given  by  (3.14)  for 
p  — *  oo  and  by  (3.15)  for  p  —  0.  The  quantity  ag  can  be  chosen  to  increase  the  sensitivity 
of  the  sensor  to  the  pressure  gradient  (a  reasonable  value  for  eg  is  5).  Crucially,  this  sensor 
will  not  switch  on  in  boundary  layers  where  the  pressure  gradients  are  weak. 


An  additional  modification  to  the  flux  evaluation  method  is  required  for  hypersonic 
flows.  In  the  stagnation  region  of  blunt  bodies,  the  convection  speeds  are  small  relative 
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to  the  sound  speed.  In  addition,  if  there  is  minor  misalignment  of  the  grid  with  the  bow 
shock  wave,  there  can  be  error  generated  by  the  bow  shock.  This  error  can  become  trapped 
in  the  stagnation  region,  resulting  in  the  “carbuncle”  phenomenon.  In  this  situation,  the 
error  overwhelms  the  actual  flow  physics  and  causes  the  bow  shock  to  lens  upstream  in 
an  aphysical  manner.  There  are  various  ways  to  prevent  the  formation  of  the  carbuncle  - 
the  best  approach  is  to  carefully  align  the  grid  with  the  bow  shock  wave.  However,  this  is 
not  always  possible,  and  an  eigenvalue  limiter  is  commonly  applied.  Here,  the  eigenvalues 
appearing  in  (3.13)  are  modified  to  prevent  them  from  going  to  zero  as  Mach  number 
approaches  zero.  For  example,  let: 

A±=  ((AUv^FjUMJ)  (3.19) 

Where  ee  is  often  taken  to  be  about  0.3.  This  eigenvalue  limiter  reduces  the  build-up 
of  error  in  the  stagnation  region  and  helps  prevent  (but  does  not  always  eliminate)  the 
carbuncle  from  forming.  It  is  important  to  make  ee  =  0  in  the  wall-normal  direction 
because  it  can  cause  artificial  diffusion  of  the  boundary  layer. 


This  approximation  to  the  flux  is  only  first-order  accurate  in  space  and  is  essentially 
worthless  for  predicting  heat  transfer  rates  for  hypersonic  aerothermodynamics  simulations. 
There  are  many  approaches  to  obtaining  higher-order  accuracy  for  conservation  laws.  In 
these  notes,  we  discuss  one  such  approach  that  has  been  shown  to  be  effective  for  a  wide 
range  of  hypersonic  flows.  It  compares  favorably  with  other  popular  upwind  approaches 
(Ref.  10),  and  its  form  is  easy  to  linearize  for  use  in  implicit  methods. 


The  key  issue  associated  with  obtaining  second-order  accuracy  is  how  to  accurately 
project  the  cell-centered  data  to  the  faces  without  introducing  numerical  problems.  In  the 
approach  above,  the  upwind  cell-centroid  value  of  U  is  taken  as  the  element  face  value  for 
use  in  the  flux  expression.  To  obtain  second-order  accuracy,  we  require  a  linear  fit  to  U 
for  a  more  accurate  value  for  U,<  i . 

2 

We  use  the  MUSCL  approach  developed  by  van  Leer  (Ref.  65)  as  our  primary  approach 
for  approximating  U  at  the  element  face.  A  simple  upwind  extrapolation  of  the  conserved 
variables  to  the  face  on  a  uniform  grid  would  result  in  a  flux  of  the  form: 


uh  =  lu<  ~  \u<-' 

U,R+i  =  \uw  -  \ui+2 


(3.20) 


and 


FU 


=  A'+ !  U 


^+2  ^  2 


(3.21) 
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Where  UL  and  UR  represent  the  approximations  to  Ui+ 1  using  left-  and  right-biased  data. 
(For  simplicity  we  have  assumed  a  uniformly  spaced  grid  indexed  by  i.) 

This  approach  yields  a  formally  second-order  accurate  flux.  However,  it  will  cause 
problems  near  strong  gradient  regions  because  the  extrapolation  can  result  in  aphysical  flow 
states.  Therefore,  these  gradients  must  be  sensed  and  the  extrapolation  reduced  or  caused 
to  revert  to  the  first-order  flux  above.  With  the  MUSCL  (monotone  upwind  schemes  for 
conservation  laws)  approach,  the  extrapolations  are  limited  to  prevent  spurious  oscillations 
near  large  gradients.  For  example,  the  extrapolation  of  each  variable  in  U  could  be  “slope 
limited.”  For  example,  consider  limiting  the  variables  (j)L  and  (j)R  with: 

fif+i  —  0*  +  x  Ihn  ((j)i+ 1  —  4>i  —  4>i-i) 

1  1  (3.22) 

—  4>i+ 1  —  2  (^*+2  _  4>i+li  0j+l  —  0i) 

where  the  limiter  function  can  take  many  forms.  We  use  a  minmod  limiter  that  takes 
the  minimum  (in  magnitude)  of  the  two  arguments  if  they  have  the  same  sign;  otherwise 
its  value  is  zero.  Note  that  this  approach  takes  the  smaller  of  the  two  possible  changes 
to  (j)  when  the  sign  of  the  slopes  is  the  same,  and  uses  no  second-order  correction  when 
the  slopes  are  of  different  sign.  Such  an  approach  can  be  shown  to  be  total  variation 
diminishing  (TVD)  for  the  linear  wave  equation. 

The  most  obvious  method  would  be  to  use  the  MUSCL  approach  on  each  of  the 
conserved  variables.  Thus,  each  conserved  variable  is  slope-limited  and  extrapolated  as 
above.  However,  we  have  found  that  a  more  robust  and  accurate  second-order  extrapolation 
method  can  be  obtained  by  applying  the  MUSCL  approach  to  the  primitive  variables  and 
then  constructing  the  conserved  variables  from  those  quantities.  With  this  approach, 
we  compute  p^’R1uL,R1vL,R1e^,R:  and  pL>R  using  the  above  expressions  and  form  UL,R 
from  these  quantities.  This  form  of  second-order  fluxes  is  recommended.  It  gives  more 
accurate  results,  and  is  significantly  more  robust  for  large  time  steps  than  the  simple 
upwind  extrapolation  (3.19)  or  the  MUSCL-based  conserved  variable  extrapolation. 

A  further  issue  involves  how  to  extrapolate  fluxes  on  non-regular  grids.  The  pre¬ 
ceding  discussion  assumes  that  the  data  are  available  from  regularly  spaced  neighboring 
elements  so  that  the  extrapolations  can  be  performed.  However,  on  a  general  triangu¬ 
lar/quadrilateral  grid  (tetrahedral/prism/pyramid/hexahedral  grid  in  three  dimensions) 
this  is  not  the  case,  and  there  is  not  a  single  neighboring  element  in  a  sensible  upwind 
direction.  In  this  situation,  it  is  necessary  to  perform  a  more  sophisticated  gradient  cal¬ 
culation  using  a  cloud  of  neighboring  points.  There  are  many  possible  ways  to  form  this 
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gradient  (Ref.  38),  but  we  have  found  that  a  weighted  least-squares  approach  gives  the 
most  accurate  results.  Here  planes  (hyperplanes  in  three-dimensions)  are  fitted  through  a 
relevant  cloud  of  nearby  data  points,  and  the  slope  of  each  variable  is  computed  from  the 
slope  of  the  plane.  This  slope  is  then  used  in  the  MUSCL  limiter  function  shown  above. 
To  be  specific,  consider  Figure  3.3  which  shows  an  example  for  which  there  is  sufficient 
information  to  form  the  gradients  with  element-centroid  data.  Thus,  the  expressions  above 
are  used  to  construct  the  fluxes  for  this  case.  Figure  3.4  shows  an  example  in  which  three 
cell-centered  values  are  available  to  construct  the  gradients  used  in  the  MUSCL  slope  lim¬ 
iter.  In  this  case,  the  weighted  least-squares  approach  is  used  to  evaluate  the  gradient  at 
the  left-side  element  centroid.  Then  the  extrapolated  variable  in  this  case  is: 

(j)L  —  fa  +  Vfa  ■  dr  (3.23) 

where  fa  is  the  variable  in  the  left  face  neighbor  element,  Vcy  is  its  gradient  there,  and  r 
is  the  face-to-centroid  vector. 


Fig.  3.3  -  Example  showing  how  element- centered  data  are  used  to  construct  the  flux  at  the 
highlighted  face. 
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Fig.  3.4  -  Example  showing  how  three  element- centered  values  and  one  weighted  least- 
squares  gradient  are  used  to  construct  the  flux  at  the  highlighted  face. 

3.2.2  Evaluation  of  the  Diffusive  Fluxes 

For  a  structured  grid  with  ordered  i.j  elements  (or  i.  j.  k  in  three  dimensions),  it 
is  possible  to  use  grid  metrics  to  evaluate  the  viscous  fluxes.  This  is  straight-forward, 
and  simply  involves  computing  unbiased  gradients  of  the  relevant  variables.  For  a  general 
unstructured  grid,  a  different  approach  must  be  used.  Two  approaches  are  commonly 
used  in  the  literature:  either  the  Green- Gauss  theorem  is  used  to  evaluate  a  gradient  by 
summing  around  the  surface  of  the  element,  or  a  weighted  least-squares  approach  is  used. 
We  favor  the  latter  approach,  though  both  are  inaccurate  in  regions  of  large  grid  stretching 
(which  is  precisely  where  accurate  gradients  are  required).  This  problem  is  particularly 
severe  in  regions  of  high  cell-aspect-ratio  (CAR),  which  is  the  ratio  of  the  longest  element 
side  to  its  shortest  side.  For  this  reason,  we  use  a  deferred  correction  approach  (Ref.  24) 
that  corrects  the  gradient  estimate  using  the  data  nearest  to  each  face. 

At  an  element  face,  the  gradient  of  some  variable  d>  can  be  written  as 

Vcj)  —  (yf>  ■  n)h  +  (V(/>  —  (V</>  •  n)n)  (3.24) 

Now,  using  the  terminology  illustrated  in  Fig.  3.5,  we  can  correct  the  face  gradient  estimate 
using  the  values  of  <fL,R 

=  (e  •  h)h  +  ^  (/  -  n  0  n)  ((V^)L  +  (V</>)fl)  (3.25) 

where  ( V<p)L,R  are  the  weighted-least  squares  gradient  estimates  at  the  left  and  right  cell 
centers.  This  approach  significantly  improves  the  gradient  values  in  high  CAR  regions. 
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Fig.  3.5  -  Deferred  correction  nomenclature. 


3.2.3  Diagonalization  of  the  Flux  Jacobian 


Now,  let  us  discuss  how  to  diagonalize  the  Jacobian  matrix,  A.  (Here  we  use  A 
for  simplicity;  the  diagonalization  of  A'  follows  trivially.)  The  straight-forward  approach 
would  be  to  form  A  and  then  find  the  eigenvalues  and  eigenvectors.  This  is  complicated 
and  difficult  to  do.  It  is  easier  to  diagonalize  A  using  a  different  set  of  variables  and  then 
transform  back  to  the  conserved  variables.  A  good  choice  is  the  vector  of  “primitive” 
variables 

T 

V  =  (pi,  p2,  •••,  pns,  u,  v,  ev,  p)  ,  (3.26) 


where  ev  —  Ev/p.  Then,  we  can  write 

_  dF  _  dU  dV  dF  dV 
~  dU  ~dV  dUdV  dU 

It  turns  out  that  it  is  easier  to  diagonalize  the  matrix  ^7  than  A  itself. 


(3.27) 


We  can  compute  these  matrices,  but  we  need  some  intermediate  results,  namely  deriva¬ 
tives  of  p  with  respect  to  the  conserved  variables,  and  derivatives  of  E  with  respect  to  the 
primitive  variables.  We  can  write  p  in  terms  of  U  as: 


and  therefore, 


P  = 


Yd s  Ps  Ma 
Yds  Pscvs 


E  —  Ev 


(3.28) 


(3.29) 
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Where  we  have  defined 


Cqj  - 


E 


-'VS  1 


*=z 


psR 

pMs 


(3.30) 


The  derivatives  of  E  with  respect  to  the  non-conserved  variables  are  computed  by 
writing  E  in  terms  of  these  variables  as 


e  —  ^ 6  R  p + ^2  Psev  +  \  ^2  Ps(u2 + t’2)  +  ^2 
l^sPsJT  as  s 


The  derivatives  that  result  are 

dE 


LSJ 

dE 


cvR 

dps  ~  y°vs  ~  MsR) 
dE 

du 
dE 

de„ 


^\T  +  ev  +  T  ^2)  T  fig: 


=  pu, 


=  Pi 


=  pv, 


dv 
dE  cv 

dp  R 


(3.31) 


(3.32) 


The  Jacobian  matrix,  that  appears  in  (3.27)  may  be  constructed  from  these 

derivatives: 

■  0  pi 


( u  0  . 

0  U  ...  0  p2 


dV  dF 
dUdV 


\ 


0  0  ...  U  pns 

u  00  1/p 

0  u  0  0 

0  0  u  0 

pa 2  0  0  a/ 


The  speed  of  sound,  a,  has  been  defined  such  that 


pa2  = 


dP  ,  dP  ,  dp  dp  ,  sdp 

E  ^  +  "sk  +  (E  +  p)  be’ 

which  may  be  simplified  using  the  derivatives  given  above  to  the  expression 

—)rt 


a  =  1  + 


Cqi  ' 


=  7i2T. 


(3.33) 


(3.34) 


(3.35) 


Where  we  have  defined  7  to  be  the  ratio  of  the  frozen  translational-rotational  specific  heats 
of  the  gas  mixture. 
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It  is  straight-forward  to  diagonalize  the  Jacobian  (3.27).  If  we  write 

dV  dF 
dUdV 

the  eigenvalues  are: 


CA  = 


(l  0 

0  1 


0  0 


V 


=  C^\ACA, 

\T 

(3.36) 

,  u,  u  +  a,  u,  u, 

s 

its 

u  —  a)  , 

(3.37) 

the  matrix  CA  is 

Pi /a2  0  0 

—ci/ a2  \ 

p2/a2  0  0 

-c2/a2 

Pns/a2  0  0 

- Cns/a 2 

(3.38) 

pa  0  0 

1 

0  1  0 

0 

0  0  1 

0 

—pa  0  0 

i  y 

Where  cs  =  ps/ p  is  the  mass  fraction  of  species  s. 


Note  that  since  the  equations  that  describe  the  reacting  flow  have  the  same  features 
as  the  perfect  gas  equations,  all  of  the  modern  upwind  flux  evaluation  methods  may  be 
used. 


3.2.4  Jacobian  of  the  Source  Vector 


To  achieve  good  convergence  rates,  it  is  necessary  to  exactly  evaluate  the  Jacobian  of 
the  source  vector,  C  =  dW/dU.  In  my  experience,  every  quantity  that  appears  in  IF  must 
be  differentiated  exactly.  In  many  cases,  neglecting  a  seemingly  small  term  can  change  the 
sign  of  some  elements  of  C ,  making  the  method  converge  slowly. 


There  are  many  different  approaches  that  may  be  taken  to  reduce  the  difficulty  of 
the  algebra.  For  example,  Gokgen  (Ref.  14)  explicitly  expresses  IF  as  a  function  of  the 
temperatures 


W(U)  =  W(U,  T(U),  TV(U )), 


(3.39) 


then  he  computes 

dW  dW  dT  dW  dTv 
~  ~dU  +  ~dTdU  +  ~&TV~W' 


(3.40) 


The  need  for  the  correct  linearization  of  the  source  term  cannot  be  understated.  Even 
small  algebra  or  coding  errors  or  simplifications  to  the  linearization  can  cause  severe  prob¬ 
lems  with  numerical  stability.  Thus,  this  part  of  the  code  must  be  rigorously  checked. 
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3.2.5  Implicit  Treatment  of  the  Boundary  Conditions 

We  have  not  discussed  the  treatment  of  the  boundary  conditions.  During  the  forma¬ 
tion  of  right-hand  side  of  (3.11),  the  appropriate  conditions  at  the  boundaries  must  be 
used.  For  example,  at  surfaces  there  must  be  no  slip  at  the  surface,  the  temperatures  must 
either  be  given  by  the  isothermal  wall  condition  or  set  by  the  adiabatic  wall  condition, 
and  the  normal  pressure  gradient  must  be  zero  or  determined  from  the  normal  momentum 
equation.  The  chemical  state  of  the  gas  at  the  wall  is  found  by  the  catalytic  efficiency  of 
the  surface. 

The  implicit  treatment  of  the  boundary  conditions  is  just  as  important,  and  is  less 
straight-forward.  Consider  a  surface  i  +  |  such  that  the  element  i  is  used  to  specify  the 
boundary  condition  for  the  flux  into  element  i  +  1.  Then  we  must  express  the  change  in 
the  solution  within  the  boundary  cell,  SUi,  in  terms  of  the  change  in  the  solution  within 
the  flow  field,  SU^ fi-  This  can  be  done  for  any  boundary  condition  if  we  construct  a 
matrix  E  such  that  SUi  —  ESUi+i.  Then  we  can  absorb  the  boundary  condition  into 
the  block-tridiagonal  solution,  and  we  can  include  the  exact  boundary  conditions  in  the 
implicit  method,  resulting  in  much  improved  convergence.  For  some  boundary  conditions, 
it  may  be  difficult  or  impossible  to  find  an  analytic  form  for  E.  In  that  case,  E  can  be 
constructed  from  numerical  derivatives  of  the  wall  flux. 

3.2.6  Implicit  Viscous  Terms 

The  evaluation  of  the  viscous  fluxes  was  discussed  in  section  3.2.2  above.  We  need  to 
linearize  the  viscous  fluxes,  Fv,  for  use  in  the  implicit  method.  We  can  write  this  flux  at 
the  unknown  time  level  as: 

F™+1  =  F ;n  +  5F™  (3.41) 

We  do  not  have  to  have  a  perfect  linearization  of  Fv .  just  the  largest  terms  need  to  be 
represented.  In  most  flows,  there  is  a  direction  in  which  the  viscous  fluxes  dominate.  For 
example,  in  high  Reynolds  number  flows  the  viscous  terms  are  large  in  the  boundary  layer, 
and  the  grid  must  be  stretched  close  to  the  surface  to  resolve  the  near- wall  gradients.  In  this 
region  the  surface-normal  viscous  fluxes  are  orders  of  magnitude  larger  than  the  streamwise 
or  spanwise  viscous  fluxes.  Therefore,  we  need  only  linearize  the  normal-direction  Fv.  This 
drastically  simplifies  the  problem,  and  it  becomes  possible  to  write  SFV  in  the  form: 

SFV  ~  Mv-^-(Nv5U)  (3.42) 

where  n  is  the  wall- normal  direction. 
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In  practice,  for  general  unstructured  grids  it  is  difficult  to  identify  the  wall-normal 
direction.  Rather,  we  simply  linearize  the  flux  due  to  gradients  of  the  face’s  nearest- 
neighbor  data. 

3.2.7  Data-Parallel  Line  Relaxation 

Let  us  return  to  the  fully  implicit  method  that  we  developed  above  in  (3.11): 

SUn  +  —  \^  An5Un  -S-  ACn5Un  =  V  Fn  ■  S  +  AWn 
V  ^  V 

sides  sides 

Substituting  the  expressions  developed  above  for  A  and  F.  results  in  a  large  linear  system 
of  equations  for  5Un.  Because  we  have  used  an  upwind  method,  the  resulting  system  of 
equations  is  diagonally  dominant,  making  it  amenable  to  solution  with  iterative  methods. 
There  are  many  such  methods  in  the  literature,  and  each  has  its  pros  and  cons.  Over 
the  past  ten  years,  we  have  been  using  the  data-parallel  line-relaxation  (DPLR)  method 
(Ref.  67)  for  the  solution  of  this  system  of  equations.  This  method  is  a  parallelizable 
variant  of  the  Gauss-Seidel  line-relaxation  method  of  MacCormack  (Ref.  33,  35),  and  is  at 
the  core  of  the  NASA  DPLR  multi-block  structured  grid  code.  This  method  is  designed  for 
use  on  parallel  computers,  and  is  ideally  suited  to  the  solution  of  wall-bounded  hypersonic 
flows.  Recently,  we  have  generalized  the  DPLR  method  to  a  certain  class  of  unstructured 
grids. 

The  DPLR  approach  takes  (3.11)  and  recognizes  that  there  is  strong  physical  coupling 
in  the  surface-normal  direction.  If  a  grid  has  been  generated  that  has  lines  of  elements 
running  out  from  the  surface  (and  preferably  through  the  bow  shock  wave  to  the  free- 
stream),  (3.11)  can  be  modified  to  reduce  its  cost  of  solution.  The  implicit  terms  due  to 
the  fluxes  that  are  transverse  to  the  wall-normal  lines  of  elements  can  be  moved  to  the 
right-hand  side  and  their  influence  included  through  a  series  of  sub-iterations.  This  results 
in  a  series  of  block-tridiagonal  solutions,  rather  than  a  full  matrix  solve.  This  method  has 
the  form: 


<H7(0)  +  ^  V  An5U^  ■  S  -  AtCn5U^  = y  Fn  ■  S  +  AtWr 


on  lines 


sides 


Then  for  k  =  1,  k 
A 
~V 


SUM  +  ^  y  AnSU^  ■  S  -  At  Cn5U^  =  -  ^  y  Fn  ■  S  +  At  Wr 
V  V  ^ 


on  lines 


v 

A 

~v 


sides 


y  An5U{k~V  •  s 


(3.43) 


off  lines 
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And  finally: 

5Un  =  SUkrnax 

One  detail  is  that  we  use  only  the  nearest-neighbor  data  in  the  linearization.  That 
is,  the  terms  due  to  the  higher-order  variable  extrapolation  are  ignored  (for  example  in 
(3.21),  we  would  not  include  the  terms  due  to  Ui- 1  and  Ui+2)-  This  simplifies  the  linear 
system,  and  (as  far  as  we  know)  does  not  diminish  the  convergence  rate  of  the  method. 

This  method  can  be  implemented  efficiently  in  parallel  with  MPI  (message  passing 
interface),  and  most,  if  not  all,  of  the  communication  costs  can  be  hidden  through  asyn¬ 
chronous  communication  protocols.  Excellent  scaling  has  been  obtained  on  a  wide  range 
of  parallel  computers. 

The  unstructured  grid  implementation  of  this  method  is  the  same  as  above,  but  re¬ 
quires  the  construction  of  surface-normal  lines  of  regular  elements.  These  elements  do  not 
have  to  be  regularly  connected,  but  are  simply  identified  as  being  in  a  wall-normal  line  of 
either  hexahedral  or  prismatic  elements.  Then  the  block-tridiagonal  solver  must  be  gen¬ 
eralized  to  allow  the  solution  on  non-regularly  connected  lines.  This  is  a  straight-forward 
generalization  of  standard  block-tridiagonal  solvers. 
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4.  Examples 

4.1  Apollo  Command  Module  Flow 

The  prediction  of  hypersonic  entry  flows  is  still  a  challenge  for  the  best  numerical 
methods.  There  are  a  number  of  sources  of  potential  problems,  and  these  include: 

•  The  main  quantity  of  interest  is  the  heat  flux,  which  is  a  gradient-based  quantity  and 
is  inherently  more  difficult  to  predict  than  the  pressure. 

•  Many  hypersonic  entry  vehicles  have  a  large  stagnation  region  with  high  heat  fluxes 
and  low  convection  speeds.  Error  can  get  trapped  in  this  region  and  accumulate, 
destroying  the  solution.  In  some  cases,  resulting  in  a  so-called  “carbuncle.” 

•  The  stagnation  region  is  bounded  by  a  very  strong  shock  wave  that  can  inject  large 
error  into  the  flow  field. 

•  The  solution  is  very  sensitive  to  grid  resolution  and  grid  alignment  with  the  shock 
wave.  All  widely  used  methods  solve  the  equations  of  gas  dynamics  in  the  grid  direc¬ 
tions.  If  the  shock  wave  is  not  aligned  with  the  grid,  the  shock  will  become  stair-cased 
across  several  cells  and  large  error  will  be  generated. 

•  Interestingly,  solution  quality  and  accuracy  may  be  adversely  affected  by  grid  refine¬ 
ment  and  stretching  in  the  stagnation  region. 

•  The  solution  is  sensitive  to  the  level  of  dissipation  used  in  the  flux  evaluation  and 
limiters.  Even  subtle  changes  in  the  flux  method  can  make  large  differences  in  the 
predicted  heat  fluxes. 

In  general,  great  care  must  be  taken  in  the  grid  generation,  a  grid  topology  with  a  patch 
in  the  stagnation  region  should  be  used,  and  the  grid  resolution  should  be  as  uniform 
as  possible  in  the  stagnation  region.  The  solution  quality  will  increase  when  the  grid 
is  aligned  with  the  bow  shock  wave.  A  certain  amount  of  skepticism  in  the  results  and 
patience  is  required  to  obtain  reliable  results  for  problems  at  high  Mach  number  and  with 
large  stagnation  regions. 

Consider  for  example  the  forebody  of  a  re-entry  vehicle  similar  to  the  Apollo  Command 
Module.  This  is  a  segment  of  a  sphere,  with  a  smaller  radius  on  its  edge.  (This  geometry 
is  a  ±20°  segment  of  a  10  m  radius  sphere,  with  a  0.5  m  radius  cylindrical  leading  edge 
added.)  Figure  4.1  shows  a  possible  surface  grid  for  this  geometry.  A  two-dimensional  grid 
for  one  surface-normal  slice  of  the  domain  has  been  rotated  about  the  symmetry  axis.  This 
results  in  pie-shaped  cells  with  a  singular  axis.  Such  a  grid  may  result  in  poor  solutions 
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because  of  the  large  variation  in  flow  resolution  in  the  subsonic  region. 

Figure  4.1  also  plots  the  computed  surface  pressure  and  heat  transfer  rate  using  this 
grid  for  the  conditions:  —  8.364  x  10-5  kg/m3,  —  219.8  K,  Voo  —  7414  m/s,  a  —  19°, 

and  Tvjaii  =  1500  K.  This  calculation  was  performed  with  the  minmod  limiter  and  the 
variables  ps,  u,  v,  w,  ev,  and  T  were  extrapolated  and  limited  to  obtain  second-order 
accuracy.  Note  that  this  is  a  particularly  difficult  condition  because  the  angle  of  attack 
produces  a  very  large  subsonic  region  where  error  can  accumulate.  This  grid  has  78  points 
in  the  radial  direction,  61  points  in  the  circumferential  direction,  and  121  points  in  the 
wall  normal  direction.  No  attempt  has  been  made  to  align  the  grid  with  the  shock  wave. 


Fig.  4.1  -  Axisymmetric  grid  on  surface  of  Apollo  Command  Module-like  shape;  surface 
pressure  (left)  and  heat  flux  (right). 


Fig.  4.2  -  Patched  grid  on  surface  of  Apollo  Command  Module-like  shape;  surface  pressure 
(left)  and  heat  flux  (right). 
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A  better  grid  generation  strategy  is  shown  in  Figure  4.2;  here  a  patch  is  added  at  the 
nose,  and  a  C-grid  is  wrapped  around  the  patch.  This  results  in  more  uniform  grid  spacing 
in  the  critical  stagnation  region.  In  addition,  this  grid  has  been  adapted  to  the  bow  shock 
so  that  the  error  produced  by  the  shock  is  reduced.  We  are  now  able  to  obtain  a  reasonable 
solution  for  this  problem.  Figure  4.2  plots  the  pressure  and  convective  heat  flux  for  the 
free-stream  conditions  given  above.  Note  that  heat  flux  no  longer  has  a  minimum  near  the 
stagnation  point,  unlike  the  result  obtained  on  the  previous  grid. 

The  main  conclusion  to  take  from  this  brief  comparison  is  that  the  solution  in  the 
stagnation  region  can  be  extremely  grid  sensitive  (at  least  using  the  methods  discussed 
here).  Great  care  must  be  taken  with  grid  generation,  and  for  problems  with  large  stag¬ 
nation  regions,  the  grid  must  be  aligned  with  the  shock  wave  to  obtain  reliable  results. 
Furthermore,  it  is  recommended  that  several  numerical  flux  methods  be  used  to  assess  the 
sensitivity  of  the  results  to  the  numerical  approach. 

The  primary  source  of  error  is  associated  with  how  the  strong  bow  shock  crosses  the 
grid.  If  there  is  perfect  alignment,  the  CFD  method  exactly  reproduces  the  shock  jump 
conditions  and  there  is  no  error.  However,  when  the  grid  is  not  aligned  with  the  shock, 
a  spurious  component  of  velocity  tangent  to  the  shock  wave  is  produced.  This  error  acts 
as  a  source  of  vorticity  at  the  shock  wave,  which  can  accumulate  in  the  stagnation  region. 
Clearly  methods  that  are  less  sensitive  to  grid  orientation  are  needed  and  new  work  in 
multi-dimensional  and  rotated  Riemann  solvers  may  help  reduce  this  dependence. 

4.2  Double-Cone  Flow 

The  double-cone  flow  held  discussed  in  Section  1.1.1  is  an  interesting  test  case  for 
evaluating  numerical  methods.  In  Figure  1.1,  the  separation  zone  that  forms  between  the 
two  conical  sections  is  shown.  The  size  of  the  separation  zone  can  be  detected  in  the  heat 
flux  and  pressure  measurements,  as  seen  in  Figures  1.2  and  1.3.  The  computed  size  of  this 
separation  zone  is  a  direct  measure  of  the  quality  of  the  numerical  solution.  Coarse  grids 
produce  a  small  separation,  and  numerical  flux  functions  with  large  dissipation  also  under¬ 
predict  the  separation  length.  Thus,  numerical  methods  can  be  evaluated  by  simulating 
the  double-cone  how  and  comparing  the  size  of  the  separation  zone. 

Figure  4.3  plots  the  separation  zone  length  as  a  function  of  the  grid  spacing  for  a 
variety  of  numerical  methods.5  This  plot  shows  that  all  of  the  methods  converge  to  the 

5  This  work  was  done  in  collaboration  with  Dr.  Marie-Claude  Druguet  of  IUSTI  —  Ecole  Polytechnique 

Universitaire  de  Marseille. 


15-46 


RTO-EN-AVT-162 


Computational  Fluid  Dynamics  for  Atmospheric  Entry 


same  separation  zone  length  as  the  grid  is  refined.  However,  the  rate  of  this  convergence 
varies  with  the  method,  and  for  a  given  grid  size  the  accuracy  of  each  method  is  different. 
Furthermore,  the  effect  of  the  slope  limiter  can  be  readily  seen,  with  the  superbee  limiter 
providing  more  accurate  results  than  the  minmod  limiter,  for  example.  Therefore,  making 
two  simulations  on  different  sized  grids  will  give  an  immediate  assessment  of  the  accuracy 
of  a  numerical  method.  It  should  be  noted  that  these  calculations  were  performed  with  a 
perfect  gas  model  without  vibrational  energy  relaxation  effects.  Thus,  the  experimentally 
measured  separation  zone  is  different  that  those  plotted  here.  See  Ref.  10  for  more  details. 

It  is  important  to  note  that  these  flows  take  a  long  time  to  evolve,  and  it  is  important 
to  carefully  monitor  the  convergence  to  steady  state.  Physically,  these  flows  take  at  least 
150  flow  times  to  converge,  where  one  flow  time  is  based  on  the  free-stream  speed  and  the 
length  of  the  geometry  (Ref.  12).  Thus,  time-like  simulations  must  be  computed  for  at 
least  this  length  of  time  before  the  solution  can  be  considered  to  be  converged. 
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Fig.  4.3  -  Size  of  the  separation  zone  versus  the  square  of  the  grid  spacing  in  the  streamwise 
direction  (Ref.  10). 
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4.3  Mach  8  Waverider  Flow 

A  third  example  involves  the  simulation  of  a  waverider  geometry  at  Mach  8  conditions 
to  compare  with  experimental  measurements  made  in  the  AEDC  Tunnel  9  facility  (Ref.  9). 
These  simulations  were  performed  on  three  grids  ranging  in  size  from  2.5  to  8.5  million 
elements.  The  outer  domain  was  designed  to  contain  the  bow  shock  wave  at  all  angles 
of  attack.  Figure  4.4  shows  some  snapshots  of  the  grid  and  waverider  geometry.  Two 
conditions  were  studied:  a  low  Reynolds  number  of  14.32  x  106/m  and  a  high  Reynolds 
number  of  53.84  x  106  /m. 

Natural  transition  occurs  on  the  waverider  at  the  lower  Reynolds  number  condition. 
We  did  not  attempt  to  model  the  transition  process,  but  rather  ran  fully  laminar  and  fully 
turbulent  flows  and  compared  those  results  with  the  experimental  data.  In  the  high  Re 
case,  transition  occurs  very  near  the  leading  edge,  and  therefore  we  ran  fully  turbulent  only. 
Here,  we  use  the  Spalart-Allmaras  RANS  model  with  the  Catris-Aupoix  compressibility 
correction  (Refs.  58,  7). 

To  a  large  extent  the  comparisons  are  very  favorable,  with  the  CFD  matching  the 
aerodynamic  data  across  the  angle  of  attack  sweep  at  both  conditions.  For  example,  Fig¬ 
ure  4.5  shows  the  lift  and  drag  coefficients  for  the  two  cases.  We  also  made  comparisons 
with  the  individual  pressure  and  heat  transfer  gauges  on  the  model.  Figure  4.6  summarizes 
one  such  comparison  for  the  low  Re  condition  on  lower  surface  of  the  model  (the  windward 
surface  at  a  >  —5°).  Note  the  excellent  agreement  with  the  data,  and  that  the  thermo¬ 
couple  near  the  leading  edge  shows  laminar  flow,  while  all  other  measurements  indicate 
turbulent  flow. 
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Fig.  4.4  -  Medium  grid  (5.1  million  elements)  used  for  the  waverider  simulations;  every 
second  point  shown  (Ref.  9). 


Fig.  4.5  -  Lift  and  drag  coefficients  for  the  waverider  at  the  low  Re  (left)  and  high  Re 
(right)  conditions. 
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Fig.  4.6  -  Heat  transfer  rate  comparisons  on  the  lower  surface  of  the  waverider  at  low  Re. 
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